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A COMPUTER PROGRAM FOR THE SIMULATION 
OF HEAT AND MOISTURE FLOW IN SOILS 

ABSTRACT 

This document describes a computer program that simulates the flow of heat and moisture in 
soils. The space-time dependence of temperature and moisture content is described by a set of 
diffusion-type partial differential equations. The simulator uses a predictor/corrector to numer- 
ically integrate them, giving wetness and temperature profiles as a function of time. The simulator 
was used to generate solutions to diffusion-type partial differential equations for which analytical 
solutions are known. These equations include both constant and variable diffusivities, and both 
flux and constant concentration boundary conditions. In all cases, the simulated and analytic 
solutions agreed to within the error bounds which were imposed on the integrator. Simulations 
of heat and moisture flow under actual field conditions were also performed. Ground truth data 
were used for the boundary conditions and soil transport properties. The qualitative agreement 
between simulated and measured profiles is an indication that the model equations are reasonably 
accurate representations of the physical processes involved. 



A COMPUTER PROGRAM FOR THE SIMULATION 
OF HEAT AND MOISTURE FLOW IN SOILS 

SECTION 1 - INTRODUCTION 

It is not generally feasible to attempt an exact simulation of heat and moisture flow in soils under 
field conditions. The number of variables involved is large, and it is difficult to include all relevant 
transport processes. However, simulations of this type are useful when used as a laboratory to 
assess the relative importance of the various factors contributing to the heat and moisture fluxes. 
They can also be used to predict the changes in heat and moisture profiles caused by the imposition 
of particular boundary conditions or by modifying the thermal and hydraulic properties of the soil- 
water system. For example, this simulator will be used as part of a study to assess the utility of 
radiometric measurements of soil microwave emissions as an indication of the moisture profile 
below the surface (Schmugge, 1978). Emission models that compute the brightness temperature 
as a function of soil moisture and temperature profiles are currently being evaluated. To test the 
emission response to varying boundary conditions (rainfall, water table height, surface heat fluxes, 
etc.), ground truth data must be measured under the appropriate conditions to use as input to 
these algorithms, which is not always possible. However, in the initial stages of model assessment, 
the output of this simulator can be used in place of the ground truth data as it correctly models 
(at least qualitatively) the response of the profiles to the boundary conditions. 

Section 2 contains a mathematical description of the program. As an example of the kind of system 
this simulator is designed to model, Section 2.1 presents the equations that describe a particular 
diffusion model of the transport properties. Section 2.2 is an overview of the predictor-corrector 
method used to numerically integrate the equations. A derivation of the equations used to perform 
the time integration is given in Section 2.3. 

Section 3 describes some test results. The simulator has been used to solve diffusion equations 
that conform to Fick’s law and for which analytical solutions are known. Section 3.1 reports on 
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these results. In Section 3.2, the results of a day-long field simulation are compared to the ground 
truth data. The resulting agreement indicates that the simulator solves the model equations cor- 
rectly. 

Section 4 is a complete description of the software, including a baseline diagram, all required input, 
all output, and the job control language (JCL) required for execution. 
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SECTION 2 - MATHEMATICAL DESCRIPTION 


2.1 DERIVATION OF EQUATIONS 

The slow movement of heat and moisture in a porous medium such as soil can be described by 
diffusion-type equations (Nielson et. al., 1972). In the classical diffusion theory, the flux (the 
amount of substance crossing a unit area per unit time) is proportional to the negative of the 
gradient of the concentration. The proportionality factor is the diffusion coefficient. The best 
known example of this kind of flow is embodied in Darcy’s Law (Hillel, 1971). The wetness 
flux is 

q g =-K(0)VOK0) — z) (2-D 

where q„ is the flux (cubic centimeters of water per square centimeter per second, cm/sec), K(0) 
is the hydraulic conductivity (cm/sec), \p(B) is the matric potential (cm), and z is the distance from 
some reference point. The term i// - z is the hydraulic head and is the potential energy of the soil 
water (matric plus gravitational) per unit weight of water. The function is called the matric 
potential and is the energy per unit weight required to overcome the capillary and adhesion forces 
that bind the water in the soil. Because work must be done to remove water from an unsaturated 
soil, i Jj is negative. The distance z is the gravitational potential energy per unit weight. 

K and i p are functions of volumetric wetness 6 (cm 3 water/cm 3 medium). In this application, it 
is assumed that soil properties change only with depth; thus the gradient is a derivative with respect 
to z. Therefore, Equation (2-1) may be written as 

%=-*m (jf-i) (2- 2a ) 

= -K(0) j|+K(0) <2-2b) 

The second line follows from the chain rule of differentiation. 

Defining a diffusion coefficient 

D<9)=K(0)^ (2-3) 
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yields, when inserted into Equation (2-2), 

% =-D(0)*j!+K(0) (2-4) 

The first term in Equation (2-4) is the diffusion contribution to the moisture flux due to a wetness 
gradient. 

There is a large body of experimental evidence indicating that thermal gradients induce moisture 
flow (Nielsen et. al., 1 972). For example, if a uniformly moisture soil sample is enclosed in a 
horizontal cylinder and is subjected to a thermal gradient, moisture flows from the warm toward 
the cool end. As field soil temperatures are always changing, an isothermal model such as Equation 
(2-4) is not complete; a theory that treats both heat and moisture flow in soils is necessary. In the 
following description, diffusion-type expressions for both heat and moisture fluxes are presented. 
The derivation closely follows the work of Philip and De Vries (1957). Contributions to heat and 
wetness fluxes that are proportional to wetness and temperature gradients are described. The 
conservation of mass and energy is then invoked to derive the partial differential equations that 
describe the variation with time of temperature and moisture profiles. 

The diffusive flux of water vapor in a porous medium is modeled as 

^=- D atm f (M)Vp (2-5) 

where q y = vapor flux density (gm/cm 2 /sec) 

D atm = molecular diffusivity of water vapor in air (cm 2 /sec) 
f = tortuosity and porosity function 
p = density of water vapor (gm/cm 3 ) 

Equation (2-5), with f = 1, describes the diffusion of water vapor in air (Eagleson, 1970). The 
factor f represents the reduced volume available for vapor diffusion due to the soil and water and 
the obstacles to the diffusing substance presented by the soil matrix. An experimentally deter- 
mined graph of f as a function of 6 can be parameterized by a linear function (Jackson et al., 1 974). 

f(0) = a(e-0) (2-6a) 
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where e = soil porosity 

a = constant less than 1 

The diffusivity D atm is a function of temperature and can be adequately modeled by the equation 
(Kimball et. al., 1976) 



where D 0 = 0.229 (cm 2 /sec) and T is the absolute temperature. 

The gradient in Equation (2-5) is to be evaluated in terms of moisture and temperature gradients 
as these are the dependent variables of the model. This can be accomplished by using the relation- 
ship between vapor density and relative humidity: 

p = p 0 h = p 0 ex P E(^g)/(RT)] (2-7a) 

where p 0 = density of saturated water vapor 
h = relative humidity 
g = gravitational acceleration constant 

R = gas constant for water vapor =4.615 X 10 6 (eigs/gm/degrees Kelvin (°K)) 

The vapor density p 0 depends on temperature and can be approximated by (Kimball et. al., 1976) 

p 0 (T) as exp (R 0 - (Ri /T)) (2-7b) 

where R 0 = 6.0035 

R, = 4975.9 (°K) 

T = temperature (°K) 

Equation (2-7) can be derived from the laws of thermodynamics. Assuming water vapor behaves 
as an ideal gas, an expression can be readily obtained relating the vapor pressure, the temperature, 
and the chemical potential of the gas. The chemical potential and the matric potential of liquid 
water are related because they both represent the free energy of the respective phases and the two 
phases are in equilibrium. The gas density is proportional to the partial pressure. 
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The gradient in Equation (2-5) can be expressed in terms of temperature and moisture gradients 
as follows: 

^P= V (p 0 h) = p 0 Vh + hVp 0 

The derivative of h with respect to T can be computed from Equation (2-7): 

/Aah i//g /g \ /9A 
\h/3T RT 2 VRT/ \9T^ 

= Z T il+ lnh(|) (It) 


The matric potential is dependent on temperature through the surface tension of water, which is 
responsible for the capillary force that binds the water to the soil matrix. Therefore, <// is pro- 
portional to surface tension a (Philip and De Vries, 1957) and 

/T\ 9 \p fl\do 

{j)W = {F)Wr ( 2 - 10 ) 


A table giving surface tension at a pressure of one atmosphere as a function of temperature (Eagle- 
son, 1 970) can be fit with the exponential 

a(T) = a 0 exp [~y(T — 273.1 6)] (2-1 1) 

where a 0 = 75.9 dynes (dyn)/cm 
7= 2.09 X 10“ 3 (1/°C) 

T = temperature, °K 


The derivative of \p with respect to T can be computed using Equations (2-10) and (2-1 1). Equation 
(2-9) therefore is 


fr-hlnh^ + ^f) 


( 2 - 12 ) 


The 0 dependence of h is, from Equation (2-7), 
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(2-13) 


dh 

d0 


( g\ 
\RT/ dd 


= h 1 n h In 




Matric potential ip typically changes by four to six orders of magnitude over the range of wetness 
values normally found in unsaturated soils. A comparison of Equations (2-12) and (2-13) shows 
that the variation of h with 0 is much larger than the variation of h with T, at least over the range 
of temperatures found in soils (273 to 310°K). Therefore, relative humidity h in Equation (2-8) 
is considered to be only a function of Q . 


Since water vapor behaves approximately like an ideal gas, its density depends primarily on pressure 
and temperature. Therefore, p 0 can be assumed to be. a function of temperature only, with no 
dependence on 0. With h depending only on 0 and p 0 depending only on T, Equation (2-8) 
becomes 


Vp = P 0 


w ♦ h 

W \9T/ 


VT 


Inserting Equation (2-14) into Equation (2-5) and using Equation (2-13) for dh/d0 yields 

>„ i 

0,vap 


where 


3 v =-D t ,, p VT-D„„..W 


D = D a(e-0)h 

T.vap atm v 7 ' 


D 


D atm «(e-0)p o gh 


6 ,vap 


RT 


(Mr) 


(2-14) 

(2-1 5a) 

(2-1 5b) 
(2- 15 c) 


This is the sought-for diffusion expression for the vapor flux. Diffusion co-efficients D T>vap an< ^ 
D fl p (respectively vailed the “thermal vapor diffusivity” and the “wetness vapor diffusivity”) 

depend on both 0 and T. 


The liquid flux can be computed from Darcy’s Law (Equation (2-1)). The gradient of p i* 1 terms of 
moisture and temperature gradients is 


yj, = Hdt V0 + — VT 
00 9T 


(2-16) 
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(2-17) 


Equations (2-10) and (2-1 1) give the derivative of \p with respect to temperature: 


d\jj_ 

3T 


-yxfj 


Thus, the liquid flux is 

O' =K ^ z - D .»V 9 -D Titi VT ( 2 - 18 a) 

where 

D 0 4iq =K (f|) (2-1 8b) 

D Tjiq = - K 7^ (2-1 8c) 


The total moisture flux, q 0 , is the sum of the vapor and liquid fluxes: 

% = q v + Q, (gm/cm 2 /sec) . (2-19) 

This can be written in a diffusion form by adding Equations (2-15) and (2-18): 

% = - d 0 .^-D t VT + KVz (2-20a) 

where 

D e = D 0,liq +D 0,vap (2-20b) 

D T sD r,iiq +D x,vap (2-20c) 

The volumetric water content, 0, is the volume of moisture per unit volume of soil. Because the 
density of water is 1 gm/cm 3 , 0 also represents the mass of water per volume soil, assuming that 
the water includes the liquid phase plus the gas phase. As 0 represents the mass and q a is the mass 
flux, they are related by the continuity equation: 

d9 3 - 

dt + V ‘ q e ~ 0 (2-21) 

This is a partial differential equation involving 6 as a function of depth and time. An analogous 

diffusion equation can be derived to describe the time dependence of the temperature profile as a 
function of the soil heat flux. 
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Fourier’s heat flow equation gives the heat flow due to a temperature gradient: 

? hT =-XVT (2-22) 

where q h T is the temperature-driven heat flux (calories/cm 2 /sec) and X is the thermal conductivity 
of the medium (cal/cm/sec/°K). 


To apply this equation to heat transfer in the soil, the effective thermal conductivity of the soil- 
water-air system must be known. A generally accepted model (De Vries, 1975) gives X as a weight- 
ed average over the thermal conductivity of each soil constituent: 


X = 


f X + 2 k.f.X. + k f (X + X ) 

w w j iii a a x a vap' 

f +Zk.f. + k f 

w • li a a 


(2-23) 


where f w , f., and f a are the volumetric fractions of the liquid, various solids, and the air, respec- 
tively. (It should be noted that f w and 6 are the same, and the porosity, e, is equal to f + f .) 

The thermal conductivities of each component are X w , X. , and X a . Factors k. represent the ratio of 
the average thermal gradient in the ith soil constituent to the average thermal gradient in water. 
They also depend on the shape and orientation of the soil grains. For spheroid-shaped particles, 


the k. factors are 
1 


k. 




-1 

4 

[— 

<< 

1 

Si 





\ W / 


(2-24) 


where g. is the shape factor and is equal to 1/2 for cylinders of infinite length, 1/3 for spheres, and 
0 for disks of infinite radius. 


The weight factor k a for air can be determined from Equation (2-24), with Xj equal to the thermal 
conductivity of dry air. The air shape factor g. in this case has no physical meaning. It is usually 
treated as a variable function of wetness that must be determined for each soil type. Near satura- 
tion, its value is usually accepted to be 1/3. 

The latent heat absorbed or emitted by the soil as the wetness changes state between the liquid and 
vapor phases can be an important cause of temperature fluctuations. This heat can be included in 
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the heat flux by assuming that the vapor flux carries with it a heat flux due to the latent heat of 
vaporization that it absorbed from the soil when it evaporated. This heat flux carried by the vapor 
phase is 


q 


h,v 



where L is the latent heat of vaporization (cal/gm) and q y is the vapor flux (Equation (2-15)). 
Both thermal and moisture gradients contribute to q y and therefore contribute to q h The mois- 
ture contribution is computed by inserting the appropriate term from Equation (2-15) into the 
above equation : 


q. =LD. VO 

^n,v(0) 0,v ap 


(2-25) 


The temperature gradient contribution from Equation (2-15) is included by increasing the apparent 
thermal conductivity of the air-filled pores, where the vapor phase primarily exists. This vapor has 
thermal conductivity X yap and carries heat flux -X yap VT according to Fourier’s Law, where VT is 
the temperature gradient in the pore. However, this heat flux can also be represented by the ther- 
mal term in Equation (2-15) with porosity factor f set equal to 1 : -D atm h (dp 0 /dT)VT. By 
equating these two expressions for the same heat flux, the vapor conductivity is found to be 


X = LD h I ^ 

vap atm \ 


(2-26) 


Therefore, the total heat flux in the soil is 


^h.T + ^h,v(e) 

= -XVT -LDfl VO 

t/,vap 

where X is given by Equation (2-23) and includes the vapor thermal conductivity. 


(2-27) 


The total thermal energy per unit volume of medium at temperature T is CT, where C is the volu- 
metric heat capacity (cal/cm 3 /°K) and T is the absolute temperature. The conservation of heat 
energy leads to an equation similar to the conservation of mass for water (Equation (2-21)); 


d(CT) 

dt 


+ V- q = 0 


(2-28) 
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The volumetric heat capacity of the soil is computed as a sum over the capacities of the constit- 
uents (De Vries, 1975); 

C = S.f.C. + f C +fC (2-29) 

lit w w a a 

Fractions f., f , and f are the volumetric contents of solid, water, and air; and C. , C , and C 

i’w a 7 i 7 w ' a 

are the heat capacities of the solid constituents, water, and air, respectively. 


Equations (2-21) and (2-28) describe the time dependence of soil wetness and temperature profiles. 
In this application, only vertical fluxes are considered; this constitutes a stratified model of the 
soil. Therefore, the gradient operator can be replaced by a derivative with respect to soil depth. 
Thus, the moisture and heat fluxes are, from Equations (2-20) and (2-27), 

«. = - D »(f)- D T©- K (2 - 30) 


/dT\ /d0\ 

%=- X {W- LD e^{Tz) 


(2-31) 


The time derivatives of the moisture and temperature profiles are, from Equations (2-21) and 


(2-28), 

dfl _ dq 0 
dt dz 

dT _ 1 / dq h\ 

dt C \ dz / 


(2-32) 

(2-33) 


BOUNDARY CONDITIONS 

To solve these equations, boundary conditions must be supplied for both wetness and temperature 
at the air/soil interface and in the bottom layer of the profile. In principle, either the fluxes q g 
and q h or the variables d and T could be specified. In the simulator both heat and moisture fluxes, 
q h and q fl , are computed at the surface. In this way the effects of the environment (i.e., rainfall, 
evapotranspiration, radiation, etc.) on the profile evolution can be modeled. At the bottom of the 
profile a mixture of flux and variable boundary conditions are used. One can specify constant wet- 
ness, a downward wetness flux equal to the hydraulic conductivity of the bottom layer, or a flux 
equal to zero. The temperature in the bottom layer is held constant. 
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When both temperature and moisture profiles are modeled the surface fluxes can be found by the 
solution of the heat balance equation 

S = R + LE + H (2-34) 


All fluxes are positive downward. S is the heat absorbed by the soil, R is the net radiation flux, 
LE is the evapotranspiration heat flux, and H is the sensible heat. This equation can be written 
with the temperature at the soil surface as the only unknown variable. After finding the solution, 
the heat flux q h at the surface is set equal to S, and the surface moisture flux q 0 is set equal to E. 


The heat flux absorbed by the soil can be evaluated by using the discrete form of equation 2-22 

S =~X, (IlZls\ (2-35) 

where X, = thermal conductivity of the surface layer 
Tj = temperature at the center of the first layer 
T g = surface temperature 
Zj = depth to center of first layer 

The net radiation R is usually divided into average short and long wavelength components (Eagle- 
son, 1970): 

R = R , + R, (2-36) 

The short wave component of R is 

R S hor« = ^ 1 “ A )[ 1 "( 1 - k ) N 3 I c ( 2 - 37 > 

where A = surface short wave albedo 

N = fractional cloud cover 

k = fraction of radiation transmitted by a completely cloud covered sky 
I c = insolation at Earth’s surface for a cloudless sky 
An empirical model for I c is (Eagleson, 1970) 

I Q sina exp (— . 1 28 n/ sina), o>0 


I = 

C 


2-38 


0 


, a<0 
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where I Q = short wave solar energy flux incident at the top of the Earth’s atmosphere (.033 
cal/cm 2 /sec) 

n = air turbidity factor (n ~ 2-5) 

a = angle between the sun and the local tangent plane computed from the following. 

sina = sin5 sin 0 + cos5 cos 0 cos r (2-39) 

5 = angle between Sun and plane of celestial equator (-23 deg < 5 < + 23 deg) 

0 = local latitude 
r = hour angle of the Sun 
= W d (t- 12) 
t = hour of the day 
W d = 7r/ 1 2 (rad/hour) 

The contribution to the net radiation from the long wavelength part of the spectrum is modeled by 


R, =oET 

long a a 


4 _ 


oT 


(2-40) 


where T a = air temperature 

E a = emissivity of the air 
T = surface temperature 

S 

a = Stefan-Boltzmann constant 
The air emissivity is modeled as (Eagleson, 1970) 

E = .74 + .005 e , (2-41) 

a a 

where e a is the vapor pressure in milli-bars. The vapor pressure and air temperature are supplied 
by the user, and should be taken from the same height above the surface. Equation 2-40 is a 
mathematical statement of the assumption that both the air and soil surface radiate with emissivi- 


ties of E a and one respectively. 

A standard model for the heat carried by the evapotranspirational flux is (Eagleson, 1970) 

7 P LU a 


LE = - - 


F/ 2 In 2 


© 


( e s-&>) 


(2-42) 
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where L = latent heat of vaporization 

P = density of air = 1.15 X 1 0 -3 gms/cm 3 
p = atmospheric pressure in units of e s and e a (~1000 millibars) 
V = Von Karman’s constant (2.5) 

Z o = surface roughness parameter (~ .025 cm. for smooth surface) 
U = wind velocity at height Z (cm/sec), averaged over 1 hour 
e a = vapor pressure at height Z 
e s = vapor pressure at soil surface 

7 = psychrometric constant (0.61 808 mb/°K) 

The surface vapor pressure is computed from 


e sat - saturation vapor pressure 


P„ t (T ) R T 

— sat s J gas s 


1000 


mb 


R gas = gas constant 

- density of saturated water vapor (Equation 2-7b) 
h = relative humidity of surface soil water (Equation 2-7a) 


(2-43) 


The evapotranspiration model outlined above has proven to be most accurate for computing fluxes 
which have been averaged over a time period on the order of an hour. However, in this simulator 
the model is used to calculate instantaneous fluxes. The comparison between simulated and mea- 
sured evaporation rates in Section 3 shows that this approach gives qualitatively correct rates. 

The mathematical model of Equation 2-42 is implemented as follows: 

LE = “ C o ~ c l U a (e s -e a ) (2-44) 

The user supplies the constants C Q , C p and U a and e a as functions of time. C Q would normally 
have the value zero. To input a constant rate offset in cm/sec, C 0 would equal the rate times the 
latent heat of vaporization (586 cal/cm 3 ). By comparing equations 2-42 and 2-44, it can be seen 
that Cj is equal to the following expression: 
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C, 


7PL 

pvJ '" 2 ( t ) 


The sensible heat flux H in Equation 2-34 is calculated by the Bowen ratio method. This ratio, 

0, is the ratio of sensible heat to the evapotranspiration flux. Under the assumption that the turbu- 
lent transfer coefficients for the two processes are equal, this ratio is 



(2-45) 


In the simulator the Bowen ratio is computed from Equation 2-45, and the sensible heat flux is 


then calculated using 


H = 0 LE 


(2-46) 


The terms of the heat balance equation are all functions of the unknown surface temperature T s 
and the known meterological variables e a , U a , and T a . The method for specifying the values of air 
vapor pressure, wind speed, and air temperature is described in Section 4. 

To solve the heat balance equation for T s , one starts by rewriting Equation 2-35 as (Hillel, 1977) 

Z, 

T s = x7 s+Ti 

Inserting the right hand side of the heat balance equation (2-34) for S gives 

T = Ti (R(T) + LE(T ) + H(T )) (2-47) 

The dependence of the flux terms on T s has been explicitly noted here. This equation is of the 
form 

T = F(T ) (2-48) 

and can be solved by the method of successive of approximations. A trial value for T. is chosen, 
F(T ) is evaluated, and a new value for T s results. This procedure can be repeated until satisfactory 
convergence is obtained. In the simulator the air temperature T a is used as a start value, and a 
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maximum of five iterations are allowed. Tests have shown that the process converges after one or 
two iterations. Convergence is defined as the absolute value of the change in T s between iterations 
being less than 0.1 °K. The process will always converge if the magnitude of the derivative of F 
is less than 1 , 


This derivative is a complicated function of changing meteorological variables, so an analytical study 
of the conditions required for convergence is not feasible. However from Equation 2-47 it is clear 
that this derivative is proportional to Z 1 , the depth to the center of the first layer. If the iterative 
procedure does not converge, then the program should be executed with a thinner surface layer. 
Lack of convergence will be evident because T will be unphysical, either too large or too small. 


Periods of rainfall can also be modeled. The user supplies the number of rain storms, the start and 
stop times of the rain (t Q and tj ), and the total accumulation (r tot ) for each one. A constant rate 
throughout each time interval is assumed and calculated as 


r = r ' ot 


t, -t 


(2-49) 


The short wave attenuation factor during each rain storm must also be supplied. This number is 
equivalent to the cloud attenuation factor 1 - (l~k) N in Equation 2-37. 

During periods of rain, the evapotranspiration and sensible heat fluxes are set equal to zero (LE 
and H, Equation 2-34), and the wetness flux at the surface q 0 , is set equal to the rain rate, Equa- 
tion 2-49. 


It is possible to remove all temperature dependence from the simulation. (See the description of 
the NAMELIST parameter ITEMPS). In this case the temperature profile, soil heat fluxes, and 
atmospheric heat fluxes are not modeled. This simulation model is useful when studying moisture 
flow in relatively moist soils. In this case the temperature gradient contribution to the moisture 
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flux is negligible, and execution time is significantly decreased when the temperature profile is 
not modeled. However, the evapotranspiration flux must still be estimated. To do this a Gaussian 
function of time is supplied, 

E(t) = - E exp [- k (t - 1 ) 2 ] (2-50) 

vy max r 1 v max 7 

where t is the time of maximum demand and E is the rate at this time. The variable k 

max max 

which determines the width of the Gaussian can be related to the integrated daily rate Eday as 
follows: 

OO 

E d.v - a E„„ Jdt exp [-k(t-t m J’ ] (2-51) 

day 

= E \ 

max V 


This gives 



(2-52) 


The user supplies t max , E max , and E day . The simulator computes k from Equation 2-52, and then 
Equation 2-51 is used to model the evapotranspiration flux. 


For some simulations it may be simpler to specify the integrated daily total and the width of the 
Gaussian. The maximum rate E max can be computed from these two. The exponential slope k 
equals 1/t 2 , where t is the time interval between the maximum rate and the time when the rate 
falls to 1/e of this value. Setting equation 2-52 equal to 1/t 2 and solving for E max gives 

1 E h 

£ = day 

max ^ ^ 

Therefore the user can compute the value of the required input E from E, and t . 

Figure 2-1 compares the model of equation 2-50 to measured data. The evaporation rates were 
measured during a field experiment performed in Phoenix, Arizona in 1971 (Jackson, R.D.). Sec- 
tion 3.2 of this document discusses this data in more detail. Curves a and b represent evaporation 
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rates three and ten days after irrigation. The squares and circles are the data, and the solid lines 
represent the model of equation 2-50. The maximum rate E and associated time t and the 

max max ’ 

daily total E day used in evaluating equation 2-50 were estimated from the data. 



Figure 2-1. Comparison of measured and modelled evaporation rates. The solid lines are the model 
of equation 2-50, and the squares and circles are the data. Curves a and b represent 
evaporation three and ten days after irrigation, respectively. 
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Figure 2-1 shows that equation 2-50 provides a qualitatively acceptable representation of evapo- 
transpiration during daylight hours, when the process is most important. The model does not in- 
clude the rise in the data from sunset to sunrise. Also the data points are not exactly symmetri- 
cal about the maximum. 


It is also possible to model the surface temperature and the heat balance equation (and thereby 
include the effect of the meteorological variables on evapotranspiration) without modeling the 
soil temperature profile. The surface temperature T s and average subsurface temperature T are 
modeled by the force restore method (Lin, 1980). The mathematical model is 


dT s _ 2S 2 tt rT 

dt a r 


dT S 
dt a yj 3657 r 


where S is the heat flux absorbed by the soil and 


(2-5 3a) 
(2-53b) 


a — 



In this expression X is the thermal conductivity of the surface layer, c is the heat capacity, and t 
is the number of seconds in a day. The thermal conductivity and heat capacity are computed using 
Equations 2-23 and 2-29. 

Since Equation 2-53 gives the time dependence of the surface temperature explicitly, T s is made 
one element of the state vector and is therefore known. Therefore, no iteration is required to 
solve the heat balance equation. The terms R, LE, and H, are evaluated using the state vector 
value for T , and S is computed from Equation 2-34. 

ROOT MODEL 

A model of soil water depletion by plant roots has been included as an extra term in the equation 
of continuity, Equation 2-32: 
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(2-54) 


d 0 _ _ dq 
dt dz 


“ Q (z,t) 


The sink term Q (1/sec) in Equation 2-54 is positive when water flows from the soil to the plant. 
The mathematical model is (Hillel, 1 977 ; Gardner, 1 964) 


<h(0,z)-4> (t) 

Q(z ’ t) = ~ % 

s P 

where 4> s (0,z) is the total potential energy of the soil water; 

<h(0,z)= <//(0)-z 

where i^(0) = matric potential 

-z = gravitational potential (z is the depth below the surface and is positive) 


(2-55) 


(2-56) 


The plant potential <f> p (cm) in Equation 2-55 varies with time but is assigned the same value 
throughout the root system. The soil resistance fl s (cm-sec) is inversely proportional to the soil 
conductivity and the amount of active roots: 

= 1/(B K (0) P (z)) (2-57) 

where B = constant 

K(0) = hydraulic conductivity (cm/sec) 

P (z) = relative root density (1/cm 2 ) 


The resistance to flow in the roots £2 p is also modeled as inversely proportional to the root density 
and root conductivity. The inverse of conductivity, called the specific resistance, is sometimes 
used. The plant resistance is 

fi p (z) = r/P (z) (2-58) 

where r = specific resistance to flow in the roots (sec/cm). 

Using Equations 2-57 and 2-58 for the resistances and Equation 2-56 for the soil water potential 
energy in Equation 2-5 5 and rearranging gives 
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BKP[i// - z - 4> p ] 
Q= 1 + BKr 


BKP [i//-z-$ p ] 
1 + ft /ft 

p ' s 


- (2-59) 


The important model parameters are the relative root density P (z) and the ratio of the resistances 
ft /ft . No loss of generality results from setting B = 1 , since its value can be absorbed into the 

P s 

definitions of P and r. Since Q is proportional to P, multiplying P at all depths by a constant 
would only change the rate at which the moisture profiles evolve. Since P has the dimensions 
(1/cm 2 ), it is commonly thought of as the length of active roots per volume of soil. As yet there 
is no experimental evidence that this is true; the model only requires that P (z) represent the rela- 
tive ability of the roots to absorb water at each depth. The plant potential <b p , commonly referred 
to as the crown potential, is modeled as a response to an atmospheric evapotranspiration demand 
function. 


The discrete model of the sink term as used in this simulator is 


.5^rrti 

J ‘ 1 + rK. 

J 


(2-60) 


Qj is the value of the sink in the j th soil layer, and z. is the depth to the center of this layer. K. 
and ipj are the hydraulic conductivity and matric potential of the soil water in the layer. The rela- 
tive root density in each layer Pj and the specific resistance of the plant roots r are input parameters. 
The crown potential $ p (t) is modeled as a response to a known transpiration demand function 
E pS (t). The crown potential is computed by requiring that the integral of the sink terms over the 
soil profile be equal to E pC . In its discrete form, this integral is 



(2-61) 


where dz^ is the thickness of the j th layer and N is the number of layers in the profile model. Using 
Equation 2-60 for Q. and solving for the crown potential gives 
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(2-61) 


N 

E P« (t > + I KjP^-^dz. 

j=l 

$ ( t ) = 

p N 

/ K.P. dz. 

Lj 111 
j=l 

Both E pE and \ p. are negative, so <h p is is also negative. Its magnitude can be large if either the 
demand is large, the soil is dry (so | is large), or both. The magnitude of 4> p must be less than 
the magnitude of the wilting point, 4 ) w , which is the largest potential for water intake the plant can 
create before wilting. Thus, the crown potential must satisfy the inequality 

(2-62) 

If 4> p < < h w the simulator will set <h p = <J> w . Once <b p is calculated the sink term can be evaluated 
for each layer using Equation 2-60. It must be positive for all layers, to correspond to flow from 
soil to roots. Any of the which are negative are set equal to zero. This procedure is used to 
accommodate experimental evidence that water flow from plant roots to the soil is negligible 
(reference 8). 

The transpiration demand E p is computed from the total evapotranspiration demand E (t). This 
is known from the solution of the heat balance equation, or from the function of Equation 2-50 
when heat fluxes are not modeled. A fraction f of the total demand will be satisfied by soil evapo- 
ration, E : 

’ s’ 

E s =fE (2-63) 

where typically f = 0.1 (Eagleson, 1970). 

The rest of the demand will be satisfied by plant transpiration E p2 ; 

E pS (t) = E -E s = (l— f)E(t) (2-64) 

This is the transpiration demand used in Equation 2-61 to compute the crown potential. 

SOIL HYDRAULIC PROPERTIES 

Both matric potential and hydraulic conductivity as functions of volumetric wetness can be 
modelled as follows for a wide range of soil types and textures, (Clapp and Homberger, 1978); 
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(2-65 a) 


K(0) = K 



2b+3 


He) = t s 



(2-65 b) 


where 0 is the volumetric wetness at saturation, K s and are the conductivity and matric poten- 

S 

tial respectively at saturation, and b is a parameter determined primarily by the soil texture. Repre- 
sentative values are 4 for sand to 1 1 for clay. This model has been implemented in the simulator. 


2.2 METHOD OF SOLUTION 

The solution space for 0(z,t) and T(z,t) consists of a variable time grid and a fixed space grid. The 
soil is divided into n layers. These layers need npt all have the same thickness. However, they 
should be small enough to enable 0 and T within two adjacent layers at a fixed time to be ade- 
quately represented by a linear function of depth. For a soil with n layers, Equations (2-32) and 
(2-33) become 2n partial differential equations: 

d0 i= _ dq 9; 
dt dz 

= K - D - D x (— ) (2-6 6b) 

q0 i 1 Vdz/ T i V dz/ 


(2-66a) 


dT i-_ 1 


dt 


/dq 


h. 


C } \ dz 
/dT : \ 






Mdz/ LE>0 - va P \dz 


(2-66c) 

(2-66d) 


If the force-restore method is used, the n equations 2-64c are replaced by the two equations 2-5 3a 
and b. 


Each equation describes the time dependence of 0 or T at a fixed depth. These equations are 
coupled, as spatial derivatives of 0 and T appear on the righthand sides. 
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An Adams-Bashforth numerical integration scheme is used to perform the time integration. The 
general equation that can be solved by this method is of the form 

3F = (2-67) 

where y is called the state vector. 


Knowledge of f(t, y) enables the generation of solutions y(t). In this application, state vector y 
contains 6 and T at various depths within the soil. Functions fare the spatial derivatives of the 

fluxes, the righthand sides of Equations (2-66), and these are known functions of z, t, 9, and T. 
The elements of y and 7 are 



i= 1, n 
i = n + 1 , 2n 


(2-68a) 


J 


_ d %_ 

dz 


i = 1, n 


2-68b) 


/-A (dq \ 

rrJ 1 “ n + 1,2n 


If the force-restore method is used, the temperature variables are replaced by the following; 

y n+ , =T s 

y n+2 = T 


(2-68c) 


The corresponding time derivatives f„ +1 and f n+2 are given by equations 2-53a and b. 

The time integration is performed via the following steps: 

1 . Estimate y(t + At), using known values y(t), f(t), ?(t - At), f(t - 2At), . . . , f^t - kAt). 
The estimate is usually called the predicted value y^ p ) (t + At) 

2. Compute f(t + At), using the elements of y< p > (t + At) from step 1 for d and T 

i i * 
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3. Recompute y(t + At), using known values y(t), f(t + A t) (step 2), f(t), . . . , f(t (k 1 )At). 
This is called the corrected value y< c) (t + At). 

In principle, steps 2 and 3 could be iterated until they both give approximately the same answer. 
However, it is more economical to reduce the integration step size At until convergence is achieved 

on the first iteration. 

The order of this method is determined by k, the number of back values of f that are used. A 
fourth order method is used in this application. The actual equations used to compute y (p) and y (c) 
are 

y( p ) (t + At)=y(t) + At [l + ±-vf(t) + V 2 f(t) + |v 3 f(t)+^ V 4 fit)] (2-69) 

and 

y( c > ( t + At) = y(t) + [25 1 fit + At) + 469 fit) + 1 09 V f(t) + 49 V 2 fit) + 1 9 V 3 fit)] (2-70) 

where functions V j fit) are linear combinations of the back values of f, defined by 

Vfit) = fit) -fit -At) 

V 2 fit) = V(Vfit)) = fit) - 2f(t - A t) - fit - 2At) 

V 3 fit) = V(V 2 f(t)) = fit) - 3 fit - At) +3 fit - 2 At) 

(2-71) 

-fit -3 At) 

V 4 f(t) = V(V 3 fit)) = fit) -4 fit -At) + 6 fit -2 At) 

-4 fit -3 At) + fit -4 At) 

These equations are derived in Section 2.3. 

To perform the time integration, the functions f., which are the spatial derivatives of the wetness 
and heat fluxes (Equations (2-68)), must be computed. 
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The soil is divided into n layers, each of thickness dz.. There are n+1 boundaries. The center of 
each layer is at depth z. below the surface. The diagram in Figure 2-1 illustrates this division. 
First, the fluxes at the n — 1 interior boundaries [q 0 , q T , i = 2, . . . , n] are calculated from the 

i i 

known values of T., and z.. The fluxes at the two boundaries [q ,q ,q ,q ] are 

U 1 n+l T 1 T n+l 

then calculated to satisfy the boundary conditions. The surface moisture and heat fluxes q and 

1 

q_r , are found from the solution of the heat balance equation. To hold the bottom temperature 
constant, the net flux into the bottom layer must be zero. The simulator does this by setting 
q T equal to q . Similarly, if constant wetness in the bottom is to be modelled, then a 0 is 

n+l n ’ M0 n+1 

set equal to q 0 . The f. are then computed using a linear finite difference representation of the 
derivative: 


dq 0 V "V 

__ i __ l+l i 

dz dz. 

1 


d % 

J 

C, dz 




(2-72a) 


(2-72b) 


The fluxes at the interior boundaries [q 0 , q T , i - 2, . . . , n] are computed as follows: 

i i 

1 . In each soil layer, compute K. (hydraulic conductivity), i//. (pressure head), (thermal 
conductivity), C. (volumetric heat capacity), D T (thermal diffusivity), and D (the coefficient 
of di^/d0 in the vapor contribution to wetness diffusivity; that is, D 0 vap = D 0 yap (d\p/d9)). The 
calculation of these quantities as functions of 6. Z t , and T. are described in Section 2.1. 

2. Compare the spatial derivative of the temperature and pressure head at the ith boundary: 



(2-73a) 



(2-73b) 
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n 


n+1 


T 



¥ 


NOTE: The center of the ith layer is at depth Z.. 

Figure 2-1 . Diagram Showing the Division of the Soil Profile Into n layers, 

Each One of Thickness dz. 
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It should be noted that in Equations (2-66) the spatial derivative of the wetness, d0/dz, always 
occurs in the product (d\p/d9) (dO/dz). This is equal to d\p/dz by the chain rule of differentiation; 
hence, the derivative of 1 p, and not 9, is computed. (This means that hysteresis is not included in 
the simulations). 


3. Compute the average value of K, X, D T , and vap at the boundaries by linear interpolation 
between the values in the two adjacent layers: 


K. dz + K , dz. 

TF — i i-l i-l l 

‘ dz. + dz , 

1 1—1 


(2-74) 


Similar expressions are used for A. , D T , and D 

1 i 


6 ,vap. * 


4. Compute the fluxes. From Equation (2-66), these are 


© 


q„ = K. - (K. + Dl ) ■ 

°j i « ® .vapj y, 

n =-T _ T fT 

^ h i i \dz /. 0 ,vap. jj 


d m) 

T i \dz/ . 

(2-75 a) 

/ d 

\dz /i 

(2-75 b) 


These fluxes are used to evaluate the righthand sides of Equations (2-72) and thereby perform the 
time integration. 

2.3 DERIVATION OF INTEGRATION EQUATIONS 

In the following development, the vector nature of the state vector is ignored to simplify the 
notation. The references provide a more detailed description. (Teddington, 1958; Booth, 1957; 
Baginski et. al., 1979) 

A general k-step formula to integrate Equation (2-67) can be written in the form 

Y n+1 = h f b k f n + l + V 1 f n + ' ' • + b o f n + l-k 1 

+ [B m-l Y n +a m-2 Y n-1 + • • • + a 0 Y n + l-m ] 

where h is the integration step size and choice of coefficients a i and bj defines the method. Back 
values of Y and f plus the current value of f may be used. A multistep method is explicit if b k = 0 


(2-76) 
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and implicit if b k ¥= 0. Implicit methods require a recursive evaluation of Y n+1 , as f n+1 appears in 
the formula. Despite this drawback, implicit methods are desirable to solve nonlinear problems 
because they can usually be designed to have smaller truncation errors and better stability. The 
Adams-Bashforth integrator uses an explicit method, known as a predictor, to calculate Y n+1 . 
f is then calculated using this estimate and an implicit method, known as a corrector, is used to 
compute a refined value, Y<*> . Calculation of derivative f n+ , typically involves far more execution 
time than manipulation of the integrator equations. However, in this formalism f is only computed 
once per integration step. Therefore, the extra accuracy derived from using the corrector equation 
costs only a small increase in execution time. 

The difference between the corrected and predicted values of Y n+1 yields, in most cases, a reliable 
estimate of the error in the integration. This error estimate is monitored to ensure that the errors 
are not too large. If necessary, the step size can be reduced by half. If the error is smaller than 
some lower limit, the step size is automatically doubled, thereby saving execution time. 

The predictor/corrector equations define Y n+1 in terms of Y n and derivatives f n+1 , f n , f n-1 , etc. 
(Only one previous value of Y is used.) However, the resulting equation is most conveniently ex- 
pressed in terms of functions called “backward differences,” which are linear combinations of f n . 


The backward difference operator V acting on discrete function g n is defined by 

Vg„=8 n -g„- 1 (2- 77a ) 

Higher powers of V acting on g n are computed by successive applications of V: 

V 2 g n = V(Vg n ) = Vg n -Vg n _, =g n -2g n _, +g n _2 (2-llb) 

v 3 g n = V(V 2 g n ) = g n - 3g n _! + 3g n ^ - g n _ 3 (2-11 c) 

V 4 g n = V(V 3 g n ) = g n - 4g n _, + 6g n _ 2 - 4g n _ 3 + g nw , (2-lld.) 

A useful recursion relation connecting various backward differences at two adjacent time steps is 

V k g n+1 = V k_1 g n+ j -V k “ 1 g n (2-1%) 


2-2 1 



Equation (2-76) written in terms of backward differences of f instead of f itself is 


= a, , Y + a, , Y , 

n + 1 1-1 n 1-2 n— 1 


+ . , . + a Y 

0 n+1— I 


+ h t c o f n + C 1 v f„ + • • • + c kV k f n + d f n+1 ] 


(2-79) 


To derive the values of coefficients a, c, and d, the backward difference operator v must first be 
related to derivative operator D = d/dt. A Taylor’s series extrapolation of a continuous function is 

Y(t+ At) = exp(At D) Y(t) 


In terms of discrete function g . this is 

°n 5 


where h is the step size. 


g n = exp(hD) g n _j 


(2-80) 


Equation (2-77) can be rewritten as 

g n =(i-v)-‘ g n _, 


(2-81) 


Comparing Equations (2-80) and (2-81) yields 

exp(hD) = (1 -V)" 1 (2-82) 

or 

hD = — ln( 1 - V) = V +^~+^y- + . . . + •• • (2-83) 

Equation (2-80) can be used to compute Y n+1 given Y n because the operator hD is known in 
terms of the backward differences from Equation (2-83): 

Y n+ i = exp(hD) Y n (2-84) 


This can be manipulated as follows: 


Y n+i =Y n +(exp(hD)-l)Y n 

-Y. + (exp(hD) -])(£§) Y 


Because DY n = f n , this yields 


(2-85) 
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Y =Y th ( eXp( !S ) -'> f- 


( 2 - 86 ) 


The numerator and denominator are expanded in powers of V From Equation (2-82), 


exp(hD) = pry 


= 1 + V +V 2 + . . . + V n + • • 


(2-87) 


Using Equation (2-83) for hD and Equation (2-87) for exp(hD) in Equation (2-86) yields 


Y ,= Y + h 

n + 1 n 


v+v 2 + V 3 + • . .+v n + . . . 
v + vi + vi+...+^-V.. 
2 3 n 


This division can be carried out to any number of terms. To fourth order, it is 

1 

2 


Y£> = Y„ +h [1 +^V + f 2 V 2 +fv 3 +^V 4 ] f n 


( 2 - 88 ) 


This is the predictor equation because it does not contain f J+( . The corrector equation can be 
derived by starting with the identity 


Y = Y + vY ., 

n + 1 n v n+1 


Manipulation as performed above yields 


Y = Y +V — Y 
Y n+1 Y n hD n+1 


= Y + h T^pr f . 
n hD n+1 


Using Equation (2-83) to expand hD gives 


Y X1 = Y + h 

n+1 n 


V + SL+...+S+.. 

v A r-» 


n + 1 


If the division is carried out to the fourth order, the following corrector equation results: 


Y('> a Y n + h [ 1 - i V - f 2 V ! ^ V 3 - 7 ^ 0 V 4 1 f t 


n+1 


(2-89) 


When an integration is performed, the derivative f n+1 is known, but the backward differences 
V k f are unknown. Therefore, the recursion relation equation (Equation (2-78)) is used to 

n+1 

express Y^ in terms of v k f„ • The resulting corrector equation is 
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(2-90) 


Y n+1 - Y n + (^§ 5 ) h f n+ i + 7^0 (469 + 109 V+ 49 V 2 + 1 9 V 3 )f n 

The integrator uses the predictor (Equation (2-88)) to estimate the state at the (n+1) step, and the 
derivative f n+1 is evaluated using this first value. This derivative and previous backward differences 
are then used to refine the estimate of Y n+J using Equation (2-90). The magnitude of the differ- 
ence between these two estimates, |Y^ - Y^ |, is a reliable estimate of the corrector error, and 
the step size can be changed at this point if necessary. This is determined by testing whether the 
difference for each state vector element lies within the error window for that layer. This window 
for the ith layer is the interval 

( w i » w h )j = (Ej/ED, Ej * ED) (2-9 1 ) 

where Ej and ED are input parameters. If the differences in all layers are less than the lower limit 
(Wj ), the step size is doubled. If the difference in one or more layers is laiger than the upper limit 
(w h ), the step size is halved. 


When the integration is complete, the backward difference table is updated via Equation (2-78) 
to prepare for the next integration step. In addition, if the step size is changed, new backward 
differences must be computed to reflect the new step size. These can be calculated in terms of the 
old differences in the following way. The value p is defined to be the ratio of the new interval 
length to the old length: 


h 


p = 


h 


old 


Halving gives p 1/2, and doubling results in p = 2. In terms of the new interval, the first back- 
ward difference is 


V P f n= f n-f„- P 


= f „-' 


-ph 


old 1 


Using Equation (2-82) to express the exponential in terms of V, this becomes 
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( 2 - 92 ) 


Vp f„ = [i - « -v) p ] f n 

- [p V v ; + p(p ~ 1 3 V P - 2) V j] . . . f„ 

The mth order backward difference is found by raising the first order to the mth power. 

V” f n = [i - (1 -v) p ] m f n 

This can easily be expanded in powers of V. Because the new backward differences are linear 
combinations of the old, changing the step size can simply be effected by a matrix multiplication. 

Halving is accomplished as follows: 
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Simulator output is desired at predetermined intervals. Because the Adams method integrates by 
using the largest interval yielding the required accuracy, an integration step frequently ends beyond 
the desired output time. As a result, interpolation is necessary; however, back values of the state 
vector are not available for this purpose. A method as accurate as integration and not requiring 
the storage of previous data points is available (Glang, 1971). This method essentially consists of 
an equation for integrating a partial step forward or backward using the existing table of backward 
differences. The derivation closely follows the derivation of the predictor equation. Let 


*n+p 

P h 


Then 


V = pP hD Y 

1 n+p e 1 n 


= Y + (e phD — 1) Y 


= Y. + (eP“>-l) *§ Y„ 


= Y n+ h 


e P hD _ 1 


hD 


Expressing this as a function of V using Equations (2-82) and (2-83) gives 


Y = Y+h 

n+p n 


(i -vr p - 1 


— ln(l -V) 


Expanding in powers of V gives 


Y = Y + h 

n+p n 


, P(P+ 1)„2 , P(P + !) (P + 2) „3 , 
P 2! 3! 


? 3 

V“ V 
— + - — + 


(2-95) 


To the fourth order in V , this is 
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Y . = Y + h 

n+p n 


(2-96) 


p 2 /2p 3 +3p 2 \ 2 

[ p+ 7 v+ ( ^ ) v 


/p 4 +4p 3 +4p 2 \ 3 

1 6p 5 + 45 p 4 + 1 1 Op 3 + 90p 2 ^ 4 

v 24 r + 

V 720 ' J 


This is the required interpolation formula. 


2-33 




SECTION 3 - TEST RESULTS 


3.1 COMPARISON OF SIMULATED TO ANALYTICAL SOLUTIONS 

This program has been tested by using it to solve diffusion-type partial differential equations for 
which analytical solutions are known. Numerical solutions of the following three equations were 


generated: 


1. Constant diffusion coefficient with constant wetness boundary condition at the surface: 



o 



(3-1) 


subject to 0 (z, 0) = 0j and 0 (0,t) = 0 Q . The solution is (Eagleson, 1970) 

0(z,t) = 0 o +(0 i -0 o )erf ^v=) 

2. Constant diffusion coefficient with constant flux at the surface: 

d0 =n ( d 2 0 
dt 0 \dz 2 

subject to 0(z,O) = 0 ; and — D Q (d0 (0,t)/dz) = f Q . The solution is (Eagleson, 1 970) 

|z| 


2f, 

e(z,t) = e i + -J 


Do 1 exp L^l )_M erfc 

T p \4D 0 t) 2 ertc { 2 ^r 


(3-2) 


(3-3) 


(3-4) 


3. Concentration-dependent diffusion coefficient with constant wetness at the surface: 

d0 _ d d0 \ 

dT " dz "d T 

(3-5) 

D(0) =— (l-£n0) 

2 

subject to 0 (z,0) = 0 and 0 (0,t) = 1 . The solution is (Philip, 1 960) 

0(z,t) = exp (— z/v/D 0 1) (3-6) 

In Equations (3-2) and (3-4), erf and erfc are the error function and complementary error function 
respectively. 
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For each simulation, the soil depth was set at 60 cm and the equations were integrated for 6 hours. 
The diffusion constant, D Q , was set equal to 0.01 cm 2 /sec. The flux at the bottom was set equal 
to 0 in all cases. This boundary condition is not compatible with the equations being integrated 
(as a semi-infinite medium is tacitly assumed in the analytical solution, no bottom boundary condi- 
tion exists). However, with this value of D Q , the water incident at the surface could not appreciably 
infiltrate to a depth of 60 cm in 6 hours; thus, the simulated and exact solutions should agree 
throughout the top layers. Results are reported to a depth of 25 cm. 

The simulator was run with an absolute error tolerance of 0.005. This is the largest allowable 
magnitude of the difference between predicted and corrected values for each iteration and is a test 
of convergence. In most cases the difference between the simulated and exact solutions did not 
exceed 0.001, indicating that the simulation converged to the correct value. The only exception is 
in the top layers when constant wetness boundary conditions are used (simulations 1 and 3). The 
exact solution gives constant wetness at z = 0, but the simulator requires all layers to have a finite 
thickness. For these two cases the top two layers were given a thickness of 0.1 cm and the thick- 
nesses of the deeper layers were gradually increased to a maximum of 2 cm. Because of this finite 
surface thickness, the exact and simulated solutions in the top cm of soil differed by about 0.01. 
Agreement at other depths was within 0.005. 

For case 2, the soil was divided into layers 2-cm thick throughout. The exact and simulated solu- 
tions were consistently within 0.001 . 

Figures 3-1 through 3-6 represent the results of the simulation. Figures 3-1, 3-3, and 3-5 give 
wetness profiles at 20 minutes, 2 hours, and six hours for cases 1, 2, and 3, respectively. Figures 
3-2, 3-4, and 3-6 show wetness as a function of time at soil depths of 1 , 11, and 25 cm for the 
three cases. 
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NOTE: 

Figure 



Depth (cm) 


Curves are labeled with the simulation time. The initial value of wetness is 0.2 at all 
depths, and the surface wetness boundary value is 0.9. 

3-1 . Wetness Profiles from the Solution of the Diffusion Equation with Constant Diffusivity 
and Constant Wetness Boundary Conditions (Equations (3-1) and (3-2)) 
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Wetness 



.9 1 - 



Depth (cm) 

NOTE: Curves are labeled with the simulation time. The initial value of wetness is 0.2 at all 
depths. The flux at the surface (f 0 ) has the value 4 X 10 -4 (cm/sec). 

Figure 3-3. Wetness Profiles From the Solution of the Diffusion Equation With Constant Diffusivity 
vity and Constant Boundary Condition (Equations (3-3) and (3-4)) 
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Depth (cm) 

NOTE: Curves are labeled with the simulation time. Initially, the profile has a wetness of 0.0 
everywhere, and the surface wetness boundary value is 1 .0. 

Figure 3-5. Wetness Profiles From the Solution of the Diffusion Equation With Variable Diffusi- 
vity and Constant Wetness Boundary Condition (Equations (3-5) and (3-6)) 
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3.2 FIELD SIMULATION 

To gain confidence in the ability of the program to solve the more complicated coupled system of 
diffusion-type partial differential equations, a day-long simulation of conditions in an actual field 
was performed. Initial conditions were taken from data that were measured in an experiment 
performed in March 1 97 1 (Jackson, 1 972). In this experiment the soil was irrigated, and at half- 
hour intervals thereafter for 1 6 days temperature and moisture profiles, evaporation rates, heat 
fluxes, and other pertinent meteorological data were measured. Thermal and hydraulic properties 
of the soil type were also known. The simulation was performed for a 24-hour period starting at 
the beginning of the fifth day after irrigation. 

Subsequent analysis of this data has shown that the flow theory described in Section 2 of this 
document can provide reasonably accurate moisture fluxes at intermediate values of wetness (Jack- 
son et. al., 1974), but the heat fluxes predicted by this theory are not very accurate (Kimball 
et al., 1976). Nevertheless, these expressions for moisture and temperature fluxes are generally 
accepted to be qualitatively correct; therefore, they may be used in a test of the ability of the 
simulator to qualitatively model heat and moisture flow in soils. 

Values for the various input parameters that were used in this simulation are described below. 

Semi-log plots of the hydraulic properties of Adelanto loam (the soil used in the experiment) were 
given in graphical form (Jackson, 1972). These are reproduced here in Figure 3-7, with the scales 
changed to reflect the different units used in the simulation. Data points were estimated from the 
curves and then used to find acceptable values for the parameters in the potential and conductivity 
models. Equations 2-65a and 2-65b. The data points are the logarithms of K and |i//|. The corre- 
sponding logarithms of the model equations can be used to solve for the model parameters. To 
use conductivity and potential data at the same time, the sum of the logarithms of K and \p was 
used; 
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K 
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Figure 3-7. Hydraulic properties of Adelanto loam. The solid curves are reproduced from 
Jackson (1972). 
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(3-7) 


log K(0) + log \\ P (0)1 = (b + 3) log (j) + log |i// s l + log K $ 

The model equations as implemented will not accommodate the sudden change in the slope of the 
potential curve near saturation. (A procedure to include this behavior is described in Clapp and 
Hornberger, 1 978). Therefore, a value of 0.375 was used for saturation wetness, and calculated 
matric potentials will not be valid above this value. A linear least squares fit of the data to Equation 
3-7 gives a value for the texture parameter b and the sum log \\p s \ + log K s which is valid for wetness 
in the range (.125, .375), where both conductivity and potential are known. Values for K s and 
d/ are then determined by examining the conductivity and potential data separately. 

This procedure results in a value of 5 .2 for b. However, the resulting fit does not match the matric 
potential data very well for wetness below 0.1 . A slight decrease in the texture parameter will 
greatly improve the fit for dry conditions, yet does not greatly alter the fit for wetness greater 
than 0.1. The following values for the model parameters were chosen to match the data: 

b = 5.0 

k s = 2.8 X 10" s cm/sec ^ ^ 

i |/ s = —42.7 cm 
0 = 0.375 

s 

The points on Figure 3-7 have been calculated with these values in Equation 2-65. The RMS 
errors of the fit to the logarithms are 

(log \\p\) = 0.37 (3-9a) 

a K (log |K|) = 0.48 (3-9b) 

These errors in the logarithms are equivalent to multiplicative errors in the hydraulic properties. 
With K o (0) and i// Q (0) defined to be the model values for conductivity and potential from Equa- 
tions 2-65 a and 2-65b, the one standard deviation ranges for the error in the fit are 10 ±CT k K 0 and 
10 ±CT ^ i j / Q . Inserting the values from Equations 3-9a and 3-9b for the RMS errors gives 
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(3-1 Oa) 
(3-1 Ob) 


.43 K o (0)<K(0)<2.3 K o (0) 

.33 \p Q (0)< i// (0) < 3.0 i p o (0) 

These errors are due only to the fit to the model equations. Errors in the measured data are un- 
known and are therefore not included. 


The soil porosity (e, Equation (2-6)) was assumed to be 0.39 because this is the largest value of 0 
for which the hydraulic properties were given. Also, the maximum values of wetness on the data 
base for any soil layer at any time is 0.38, which occurred in the top layers just after irrigation. 

The volume fractions of the solid constituents of Adelanto Loam are .373 for quartz and .627 for 

clay (Kimball, et al., 1 976). The volume fractions of quartz and clay for the soil-air-water system 
are therefore 

f q = .373 (l-e) = 0.228 

f c =.627 (l-e) = 0.384 (3_1 


The thermal conductivities of quartz and clay are 21 and 7 (mcal/cm/sec/°C) respectively. Equation 
(2-24) gives, for the weighting factors, 


k q = -173 
k c = .422 


(3-12) 


The value used for the shape factor g. is 0.333. 


Equation (2-23) contains contributions to the effective thermal conductivity from all solid consti- 
tuents. In the simulator, these factors are replaced by one effective solid term with volume fraction 
f = l~e. This is accomplished by defining the effective thermal conductivity of the solids X and 
weighting factor k by solving the equations 

k(l-e) \= SkjfjAj (3-13) 

k(l-e)=2k.f. (3-! 4) 
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(3- 15a) 


Using previously stated values k., f., and these equations can be solved, giving 

k = 0.329 

T= 9.72 X 10" 3 (mcal/cm/sec/°C) (3-15b) 

The volumetric heat capacities of quartz and clay have the same value, 0.5 (cal/cm 3 /°C), which was 
used for C s in Equation (2-29) to compute the volumetric heat capacity of the soil-water system. 

The soil moisture in the bottom layer was approximately constant over a 24-hour period; con- 
sequently, the- option to hold this value constant in the simulation was chosen. (JBOT = 1 in 
NAMELIST, see Section 4.) 

Values for parameters needed for the solution of the heat balance equation were either measured or 
could be deduced from the measurements. Air temperature, vapor pressure, and wind speed were 
measured at half-hour intervals, and a linear interpolation scheme was used to provide values at 
intermediate times. The subroutine which did this is listed in Section 4. 

Values for parameters in the solar radiation flux model (Equation 2-37) are as follows: X = 33.5°, 

5 = 0.0°, A = 0.3, N = 0.0, and n = 2.0. The first two variables in this list follow from the location 
of the experiment (Phoenix) and the time of year (March 8, near the vernal equinox). Both in- 
coming and outgoing short wave fluxes were measured, and the average value of their ratio (the 
albedo) over the day was 0.3. Figure 3-8 shows both modeled and measured fluxes. The qualita- 
tive agreement is very good. The major error source is the phase difference of one-half to one hour. 

Evaporation rates were measured during the experiment. Therefore the constant Cj in Equation 
2-44 was chosen to make the model calculations match the data as accurately as possible when the 
rate was the largest. The simulated and measured evaporation rates are compared in Figure 3-9, 
with the value of 0.5 X 10 -6 (cal/cm 3 /mb) used for Cj . In the nominal simulation C Q was set 
equal to zero. The simulation was also run with the value of C 0 chosen to represent the reported 
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Figure 3-8. Simulated and measured incoming and outgoing short wave radiation fluxes for March 
8 simulation. Solid lines represent the data, and the points represent the model calcu- 
lations. 
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Figure 3-9. Comparison between simulated and measured evaporation rates. 
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uncertainty of .02 mm/hr., the resulting curve is also shown in Figure 3-9. The qualitative agree- 
ment is good, indicating that Equation 2-44 can be used to model physically realistic instantaneous 
evaporation rates. 

The results of a 24-hour simulation are presented in Figures 3-1 0 through 3-15. The qualitative 
similarity between the measured and simulated variables is evidence that the program has solved 
the model equations correctly. As previously noted, the solutions of the model equations are not 
expected to agree exactly with the data because of uncertainties in soil properties. 

Figure 3-10 shows simulated and measured wetness in the top Vi cm of soil as a function of time. 
The greatest daily variation occurs in this layer, so a comparison of the simulated and measured 
wetness there provides the most stringent test of the program. Curve a is the data, and curve b 
is the simulation using nominal values of the input parameters. The quantitative agreement between 
these two curves is poor. However, some of the discrepancy can be accounted for by uncertainties 
in the input parameters. For example, the error in the evaporation rate is conservatively estimated 
to be .02 millimeters per hour (Jackson, 1972). Curve c results from increasing the magnitude of 
the evaporation rate in the simulation by this amount, and only a small change is produced. The 
values of the hydraulic parameters can also be in error by considerable amounts (Equation 3-10). 
Decreasing the magnitudes of K and \p by the factors of .33 and .43 respectively (and therefore 
reducing the value of the wetness flux everywhere) results in curve d. The simulated wetness is 
very sensitive to changes in these variables. It should also be noted that the errors quoted include 
only those due to curve fitting and extracting numbers from the semilog plot; the errors in mea- 
suring K(0) and a re riot known and have therefore not been included. Finally, curve e shows 
the effect of changing both the evaporation rate and hydraulic parameters. In this case changing 
the evaporation rate has a greater effect on the moisture content . 

Figure 3-1 1 shows the simulated and measured surface temperatures, and Figure 3-12 shows the 
corresponding heat fluxes at a depth of 5 cm. The qualitative agreement in both cases is good. 
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Figure 3-1 0. Wetness in top 14 cm. of soil as a function of time. Curve a is the data, and curves 
b-e represent simulations with various uncertainties in the input variables included. 
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Figure 3-12. Simulated (dashed line) and measured (solid line) heat fluxes at 5 cm. depth. Fluxes 
are positive away from the surface. 



The maximum error in the surface temperature is 3°C. Errors in the heat flux can be reduced by 
calibration of the thermal conductivity model (Equation 2-23). (Kimball, et al., 1976a and b.) 
This is done either by making the air shape factors (Equation 2-24) functions of soil moisture, or 
by introducting a constant factor in Equation 2-23. However, in this simulation the value 1/3 was 
used for g a for all values of wetness, and no multiplying factor was used. 

Figures 3-13 and 3-14 show wetness profiles for three different times of the day. The simulations 
represented in figure 3-1 3 were performed with nominal values for the input parameters (curve b, 
figure 3-10); figure 3-14 shows the changes when errors in the hydraulic parameters are included 
(curve d, figure 3-10). Figure 3-15 shows temperature profiles from the nominal simulation. The 
simulations were performed to a depth of 100 cm' but little variation was found below 25 cm for 
wetness and temperature in both the simulations and the data; therefore, profiles below this depth 
are not shown. The qualitative agreement is good for moisture profiles below 7 cm. The tempera- 
ture profiles agree within 3° at all depths. 

Uncertainties in the hydraulic parameters are responsible for most of the surface moisture differ- 
ence between the simulations and the data; errors in the simulation of evaporation rate are of 
secondary importance. Surface soil moisture variations are determined by the difference in the 
evaporation rate and the recharge flux into the surface layer from below. Table 2-1 shows values 
for the integrals of these fluxes over one day for the four simulations discussed. The integrated 
evaporation rate from the data is also presented. 

The differences between columns b and d show that changing the hydraulic parameters has very 
little effect on the evaporation rate but does alter the recharge flux. Therefore large variations 
in surface soil moisture can be created by changing the hydraulic parameters because the input to 
the layer changes, but output is relatively constant. (This decrease in recharge is reflected in the 
smaller change over time of the simulated moisture profiles in Figure 3-14 when compared to 
Figure 3-13). However, if only the evaporation rate is changed (column c), both surface and 
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Figure 3-13. Wetness profiles at three times for the sample simulations. Solid lines are data, and 
dashed lines are simulations. 
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Data Source 



a 

b 

c 

d 

e 

Evaporation 

.181 

.125 

.167 

.126 

.168 

Recharge 


.157 

.196 

.118 

.147 


Table 2-1 . Daily Integrated Evaporation and surface recharge fluxes, in cm of water. 
Data 

Simulation with nominal input values 

Simulation with evaporation rate increased by uncertainty 

Simulation with hydraulic parameters reduced 

Simulation with evaporation rate increased and hydraulic parameters reduced 




recharge fluxes change by almost the same amount. Thus the effect on surface soil moisture is 
minimal. The increased evaporation in this case is supplied by moisture in layers deeper than the 

surface. 
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SECTION 4 - USER’S GUIDE 


Figure 4-1 shows the JCL needed to execute the simulator on the SACC (Science and Applications 
Computing Center, GSFC) S/360-91 and S/360-75 computers. Subroutines which perform the 
integration and create printer plots are in object module format as members of the data set 
ZB2PC.UTIL.OBJ. All other subroutines exist as a non-executable load module named ZB2PCSSM 
on the system load module data set, SYS2.LOADLIB. The linkage editor is used to combine 
these and create an executable load module. In this example NAMELIST input is entered instream 
in the GO.DATA5 DD statement. The DD card labeled GO.FTlOFOOl is for printer plots (see the 
description of NAMELIST parameter IPR). The GO.FT12F001 DD statement points to the disk 
data set which will receive the temperature and moisture profiles (see the description of NAMELIST 
parameters NDISK and IUDISK). 

Figure 4-2 is an example of the JCL used to relink a subroutine created by the user and then exe- 
cute the simulator. The first job step is a compilation of the user subroutine using the standard 
FORTRANH procedure, and the second step is the LINKGO procedure from Figure 4-1. Note 
that in the present example the NAMELIST is read from the data set Z B2PC .MARCH 1 7 .NL.DATA. 

The variables in the simulator NAMELIST are described in Table 4-1 . For each variable the type 
(real (R) or integer (I)) length (4 or 8 bytes), and default value are given. If the variable name and 
value does not appear in the NAMELIST input data, then the variable will have the value given in 
Table 4-1. Variable descriptions are also given. References to defining equations are provided 
where relevant. 

In addition to values for the NAMELIST variables, the user can provide values for the air tempera- 
ture, air vapor pressure, and wind speed as functions of time. This is done by coding a subroutine 
named METEOR. The routine which is used in the simulator is given in Figure 4-3. The calling 
sequence is as follows : 
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//* SOIL SIMULATE 

//[■INK f5|f : TD LI ^ Go iREGION.GO=300K,OUT = 8 
//LINK. SYSLI B DO DISP=SHR , DSN=ZB2PC . UTI L. OBJ 
//LINK. OBJECT DO * 


INCLUDE LOADLIB(ZB2PCSSM) 
ENTRY MAIN 


/♦ 

//GO.DATA5 DO * 
SINPUT 


NL=22,DZ=1 . 0,2.0, 3. 0,4.0, 18*5.0, WAT ER*27*0. 1 , TEMPS«27*290 . 0 , 
SATW=27*0.4 ,SATP=27*-10.0,SATK«27*1 . OE-4 , EB=27* 5 . 0 , 

TST0P=8.64D5, DTOUT = 10800.0, 

WATERR=0.002,TEMERR=2.0,ED=2.0D0,HMAX=1800.0D0, 

I ROOTS = 0 , 1 T EMPS = 1 , NWATRS=2 , NWFLUX=2 , NTEMPS=2 , NTFLUX=2 , 

SPRES= 1 .OD6,ROOTS=27*0. 1 , 

I NDXW= 1 ,2,1 NDXWF= 1 ,2,INDXT=1 ,2,INDXTF=1 ,2, 

NHCUMS= 1 , IXHCUM=1 ,NWCUMS=2 , IXWCUM=1 ,2, 

NDISK=1 ,IUDISK=12,- 

RNSTRT=8.64E4,RNST0P=1 . 728E5 , RNTOT= 0 . 0 , 

NFUNCT = 8 , 

SEND 

//GO.FTIOFOOI DD SYSOUT = 8 , DCB= ( REC FM = VBA , LRECL= 1 37 , BLKSI ZE = 7265) 
//GO.FT12F001 DD DISP=SHR , DSN*ZB2PC . PROFI LE . NORAIN . DATA ( DR YNORTS) 


Figure 4.1 . Execution JCL 


//COMP EXEC FORT RANH, PARM*XREF , OU T*8 

//SYSIN DD DISP=SHR,DSN=ZB2PC. METEOR. FORT 

//LGO EXEC LINKGQ. REGION. GO=300K ,0UT=8 

/ /LINK. SYSLI B DD DISP-SHR,DSN=ZB2PC . UTI L. OBJ 
//LINK. OBJECT DD * 

INCLUDE L0ADL1B(ZB2PCSSM) 

ENTRY MAIN 
/* 

/ /GO.DATA5 DD DI SP=SHR , DSN=ZB2PC .M ARCH1 7 . NL . DATA 

//GO.FTIOFOOI DD SYSOUT *8 , DCB* ( REC FM*VBA, LR£CL* 1 37 , BLKSI ZE * 7265) 


Figure 4.2. Relink JCL 
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TABLE 4-1 NAMELIST INPUT VARIABLES 


All subscripted variables are arrays with one element per soil layer. Subscripts run from 1 to NL 
(number of soil layers) unless otherwise indicated. The NAMELIST name is INPUT. 


Integrator Control Variables 

Name Type 

TSTOP R*8 

Default 

8.64D4 

Description 

Stop time for integration (seconds 
from start) 

IFORCE 

1*4 

1 

Force integration step size to 
remain less than HMAX (0 = no, 
1 = yes) 

HMAX 

R*8 

1.8D3 

Maximum step size (seconds) 
(not applicable if IFORCE = 0) 

H 

R*8 

1.0D0 

Initial step size (seconds) (not 
applicable if IFORCE = 1 ; in this 
case H is automatically set equal 
to HMAX/512) 

WATERR 

R*4 

1.0E-3 

Error tolerance parameter for 
wetness in soil layers (Equation 
(2-59)) 

TEMERR 

R*4 

1.0 

Error tolerance parameter for 
temperature in soil layers 
(Equation (2-59)) 

ED 

R*8 

5.0D0 

Error window parameter (Equa- 
tion (2-59)) 

IROOTS 

1*4 

0 

Include water uptake by plant 
roots in the simulation (0 = no, 
1 = yes) 

ITEMPS 

1*4 

0 

Temperature model indicator: 

0 = no temperature in model, 

1 = model soil temperature pro- 
file, 2 = use force-restore equa- 
tions to model surface and average 
subsurface temperatures 
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Name 

JBOT 

Type 

1*4 

Default 

1 

Description 

Bottom wetness boundary condi- 
tion indicator: 0 = flux is zero, 

1 = wetness is constant, 2 = flux 
is hydraulic conductivity of 
bottom layer 

Output Control Variables 
Name 
DTOUT 

Type 

R*4 

Default 

1800.0 

Description 

Output period (seconds) 

NFUNCT 

1*4 

0 

Number of wetness and tempera- 
ture profiles per plot page 
(0 < NFUNCT < 10) 

NWATRS 

1*4 

0 

Number of soil layers for which 
wetness is to be plotted as a 
function of time 
(0 < NWATRS < 10) 

INDXW 

(1,1=1,10) 

1*4 

10*0 

Indices of soil layers for wetness 
versus time plots 

NWFLUX 

1*4 

0 

Number of soil boundaries for 
which wetness flux is to be 
plotted as a function of time 
(0 < NWFLUX < 10) 

INDXWF 
(1,1= 1, 10) 

1*4 

10*0 

Indices of soil boundaries for 
wetness flux versus time plots 

NTEMPS 

1*4 

0 

Number of soil layers for which 
temperature is to be plotted as a 
function of time 
(0 < NTEMPS < 10) 

INDXT 

(1,1=1,10) 

1*4 

10*0 

Indices of soil layers for tempera- 
ture versus time plots 

NTFLUX 

1*4 

0 

Number of soil boundaries for 
which heat flux is to be plotted 
as a function of time 
(0 < NTFLUX < 10) 

INDXTF 

(1, 1=1,10) 

1*4 

10*0 

Indices of soil layers for heat flux 
versus time plots 

ITABLE 

1*4 

1 

Print tables of important variables 
(0 = no, 1 = yes) 
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Name 

IPR 

Type 

1*4 

Default 

10 

Description 

Output unit number for printer 
plots (tables are printed on unit 6) 

WL 1 

R*4 

0.0 

Lower limit for wetness on plots 
(cm 3 /cm 3 ) 

WH 1 

R*4 

0.0 

Upper limit for wetness on plots 
(cm 3 /cm 3 ) 

WFL 1 

R*4 

0.0 

Lower limit for wetness flux on 
plots (cm/ sec) 

WFH 1 

R*4 

0.0 

Upper limit for wetness flux on 
plots (cm/sec) 

TL 1 

R*4 

0.0 

Lower limit for temperature on 
plots (°K) 

TH 1 

R*4 

0.0 

Upper limit for temperature on 
plots (°K) 

TFL 1 

R*4 

0.0 

Lower limit for heat flux on 
plots (cal/cm 2 /sec) 

TFH 1 

R*4 

0.0 

Upper limit for heat flux on 
plots (cal/ cm 2 /sec) 

NWCUMS 

1*4 

0 

Number of layer boundaries for 
which cummulative wetness 
fluxes are to be computed 
(0 < NWCUMS < 10) 

IXWCUM 

(1,1= 1, NWCUMS) 

1*4 

0 

Indices of boundaries for which 
cumulative wetness is to be com- 
puted 

NHCUMS 

1*4 

0 

Number of layer boundaries for 
which cumulative heat fluxes are 
to be computed (0 < NHCUMS 
< 10) 

IXHCUM 

(1,1= 1, NHCUMS) 

1*4 

0 

Indices of boundaries for which 
cumulative heat flow is to be 


computed 


1 If the lower and upper limits for any variable are equal, the actual limits used are determined 
from the data being plotted. 
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Name 

NDISK 

Type 

1*4 

Default 

0 

Description 

Flag to indicate if wetness and 
temperature profiles are to be 
written to disk or tape (0 = no, 
1 = yes) 

IUDISK 

1*4 

12 

Unit number for output device 
(disk or tape) 

Variables Defining Soil Profile and Properties 
Name Type 

NL 1*4 

Default 

200 

Description 
Number of soil layers 
(2 < NL < 200) 

DZ(I) 

R*4 

200*1.0 

Thickness of soil layers (cm) 

WATER(I) 

R*4 

200*0.25 

Initial volumetric wetness of 
soil layers (cm 3 /cm 3 ) 

TEMPS(I) 

R*4 

200*293.0 

Initial temperature of soil 
layers (°K) 

SATW(I) 

R*4 

200*0.3 

Saturation volumetric wetness of 
soil layers (cm 3 /cm 3 ) ( d s , Equa- 
tion 2-65) 

SATK(I) 

R*4 

200*1.0E-4 

Saturation hydraulic conductivity 
of soil layers (cm/ sec) (ks, Equa- 
tion (2-65 a)) 

SATP(I) 

R*4 

200 *-10.0 

Saturation matric potential of 
soil layers (cm) (t// s , Equation 
(2-65 b)) 

EB(I) 

R*4 

200*5.0 

Texture parameter of soil layers b, 
(Equation (2-65)) 

PORSTY(I) 

R*4 

200*0.45 

Porosity of soil layers (cm 3 /cm 3 ) 

TCONDSd) 

R*4 

200*2.5E-3 

Thermal conductivity of solid 
matter in soil layers (cal/cm/ 
sec/°K) 

VHCAPSa) 

R*4 

200*0.5 

Volumetric heat capacity of soil 
layers (cal/cm 3 /°K) 

FACTKS 

R*4 

0.75 

Shape factor for soil grains (k } in 
Equation (2-23)) 
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Name 

TCONDW 

Type 

R*4 

Default 
1 .3E-3 

Description 

Thermal conductivity of water 
(cal/cm/sec/°K) 

TCONDA 

R*4 

5.967E-5 

Thermal conductivity of air 
(cal/cm/sec/°K) 

VHCAPW 

R*4 

1.0 

Volumetric heat capacity of water 
(cal/cm 3 1 ° K) 

VHCAPA 

R*4 

3.0E-4 

Volumetric heat capacity of air 
(cal/cm 3 /°K) 

FACTKA 

R*4 

1.4 

Air weighting factor in definition 
of thermal conductivity (k a in 
Equation (2-28)) 

ALPHA 

R*4 

0.667 

Tortuosity factor (a, Equation 
( 2 - 6 )) 

GAMMAO 

R *8 

2.09D-3 

Surface tension temperature 
coefficient ( 7 , Equation (2-1 1)) 

RHOVPO 

R*4 

6.0035 

Constant in exponential to com- 
pute the density of water vapor 
(Equation (2-7b)) 

RHOVPT 

R*4 

4975.9 

Coefficient of 1/T in definition of 
the density of water vapor (Equa- 
tion (2-7b)) 

datmo 

R*4 

0.229 

Diffusion coefficient of water 
vapor in air at 0° C (cm 2 /sec) 
(Equation (2-6b)) 

LHEAT 

R *8 

586.0D0 

Latent heat of vaporization of 
water (cal/gm) 

THMIN 

R *8 

0.05D0 

Minimum value of soil moisture 
to support evaporation. If the 
surface soil moisture is below this 
value, evaporation is limited to 
flux into the surface layer from 
below 

SFRAC 

R *8 

0.1D0 

Fraction of evapotranspiration 
demand which is satisfied by 
evaporation (f. Equation (2-63)) 

•ROOTS(I) 

R *8 

200*0.0D0 

Root density profile (P j; Equation 
(2-60)) 
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Name 

CPMIN 

Type 

R*8 

Default 

-1.5D4 

Description 

Limiting value of crown potential, 
cm (0 w , Equation (2-62)) 

SPRES 

R*8 

1.0 D6 

Specific resistance of roots, (sec/ 
cm) (r, Equation (2-58)) 

Variables Defining the Processes at the Air/Soil Interface 


Name 

EMAX 

Type 

R*8 

Default 

3.0D-5 

Description 

Maximum rate for time dependent 
evapotranspiration model (cm/sec) 
(E max , Equation (2-50)) 

EMAXT 

R*8 

4.68D4 

Time of maximum rate, seconds 
since start of simulation (t 

x rn 3X 5 

Equation (2-50)) 

EDAY 

R*4 

1.0 

Total daily evapotranspiration 
(cm) (E day , Equation (2-51)) 

NRAINS 

1*4 

0 

Number of rain storms 

RNSTRT (10) 

R*4 

0.0 

Start time of rainfall, seconds 
since start of simulation (t , 
Equation (2-49)) 

RNSTOP ( 1 0) 

R*4 

0.0 

Stop time of rainfall, seconds 
since start of simulation (t, , Equa- 
tion (2-49)) 

RNTOT (10) 

R*4 

0.0 

Total rainfall accumulation, (cm) 
(r to{ , Equation (2-49)) 

ATTEN (10) 

R*4 

0.5 

Short wave attenuation during 
rainfall 

CEVAPO 

R*4 

0.0 

Constant term in evapotranspi- 
ration model (cal/cm 2 /sec) (C , 
Equation (2-44)) 

CEVAPI 

R*4 

0.0 

Coefficient of variable term in 
evapotranspiration model (Cj , 
Equation (2-44)) 

XLAT 

R*4 

45.0 

Latitude of simulation, degrees 
(0, Equation (2-39)) 

SUNDEC 

R*4 

0.0 

Declination of Sun, degrees 
(5, Equation (2-39)) 
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Name 

Type 

Default 

Description 

ALB 

R*4 

0.3 

Surface albedo to shortwave 
radiation (A, Equation (2-37)) 

CTRANS 

R*4 

0.2 

Cloud transmittivity (k. Equation 
(2-37)) 

CLOUDS 

R*4 

0.0 

Fractional cloud cover (N, Equa- 
tion (2-37)) 

TURB 

R*4 

2.0 

Air turbidity factor, (n, Equation 
(2-38)) 
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SUBROUTINE METEOR( T , TAIR , VAP20 . WIND) 
REAL*8 T 

DIMENSION TEMPS(50) ,VAPORS(50) ,WINDS(50) 
DATA TEMPS/ 

0.28036E 03, 

0.27846E 03, 

0.27796E 03, 

0.28106E 03, 

0.29806E 03 , 

0.30566E 03, 

0.30126E 03, 


0.28656E 03, 
DATA VAPORS/ 

* 0.49000E 01 , 
0.25000E 01 , 


0.77000E 01 , 
0.53000E 01 , 
0.53000E 01 , 
0.52000E 01 , 
0.46000E 01 , 


0.28056E 

03, 

0.28026E 

03, 

0.28046E 

03, 0.27976E 

03 

0.28166E 

03, 

0.27586E 

03, 

0.27806E 

03, 0.27746E 

03 

0.27896E 

03, 

0.27866E 

03, 

0. 27856E 

03, 0.27806E 

03 

0.28306E 

03, 

0.29236E 

03, 

0. 30076E 

03, 0.29846E 

03, 

0.30126E 

03, 

0.30206E 

03, 

0.30196E 

03, 0.30396E 

03, 

0.30576E 

03, 

0.30546E 

03, 

0. 30296E 

03, 0.30096E 

03, 

0.30186E 

03, 

0.30266E 

03, 

0. 30216E 

03, 0.30296E 

03, 

0.29936E 

03, 

0.29606E 

03, 

0.28896E 

03, 0.28786E 

03, 

0.28666E 

03, 

0. 28.556E 

03, 

0.28706E 

03, 0.28796E 

03, 

0.28736E 

03, 

0.28206E 

03, 

283.26E0 , 

, 283.06E0/ 


0.44000E 

01 , 

0.32000E 

01 , 

0.96000E 

01, 0.57000E 

01 , 

0.63000E 

01 , 

0.49000E 

01 , 

0.50000E 

01, 0.55000E 

01 , 

0.55000E 

01 , 

0.52000E 

01 , 

0. 55000E 

01, 0.55000E 

Oil 

0.61000E 

01 , 

0.67000E 

01 , 

0. 76000E 

01, 0.13500E 

02, 

0.59000E 

01 , 

0.56000E 

01 , 

0.57000E 

01, 0.55000E 

01 , 

0.57000E 

01 , 

0.55000E 

01 , 

0. 57000E 

01 , 0.57000E 

01 i 

0.51000E 

01 , 

0.50000E 

01 , 

0.48000E 

01, 0.47000E 

01 i 

0.44000E 

01 , 

0.44000E 

01 , 

0.59000E 

01, 0.59000E 

01 , 

0.47000E 

01 , 

0.51000E 

01 , 

0.46000E 

01, 0.43000E 

01 , 

0.45000E 

01 , 

0.51000E 

01 , 

4.8E0.4.4E0/ 



0.45000E 01 , 

DATA WINDS/ 

59.0, 68.6, 85.4, 94.8, 104.8, 115.8, 136.4, 114.8, 

111 A no n cn n ^ - * 


* 112.4, 

93.0, 

60.2, 

91 .4, 

98.2, 

* 99.0, 

125.0, 

107.2, 

104.8, 

1 12.6, 

* 40.6, 

67.2, 

90.4, 

112.8, 

139.8, 

*. 93.4, 

77.0, 

77.2, 

83.8, 

89.2, 

* 58.6, 

58.6, 

60.0, 

59.2. 

64.0/ 

RINTRP(A 

» B ) = A 

+ (B 

- A ) *DTRE L 

TREL 

DM0D(T 

, 8 . 64D4 )/ 1 800 . 

. 0 DO 

ITRL 

INT(TREL) 



10 

ITRL + 

1 



11 

10 +■ 1 




DTREL = 

TREL - 

ITRL 




74.6, 

08.4, 

162.4, 

78.0, 


88 . 0 , 

54.2, 

146.8, 

75.4, 


100 . 0 , 

60.4, 

127.4, 

68 . 0 , 


98.8, 

79.2, 
50.6, 

102 . 0 , 

62.2, 


TAIR 
VAPZO 
WIND 
RETURN 
END 


RINTRP(TEMPS( 10) , TEMPS ( II ) ) 
RINTRP( VAPORS ( 10) , VAPORS ( II ) ) 
RINTRP(WINDS( 10) , WINDS ( II)) 


Figure 4-3. Example of subroutine METEOR 
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T - input - Time in seconds since start of simulation (R*8). 

TAIR - output - Air temperature at reference height above the surface, degrees K, (R*4). 

VAPZO — output — Air vapor pressure at reference height, millibars (R*4). 

WIND - output - Wind speed at reference height, cm/sec (R*4) 

The output variables are used to solve the heat balance equation. In this example data are provided 
at half-hour intervals, and a linear interpolation sceme is used to provide values at other times. The 
DMOD function causes the simulator to reuse the data in case the simulation extends beyond 24 
hours. 

As long as the calling sequence convention is followed and the data have the correct units, the user 
is free to use any scheme to generate values. The relink JCL of Figure 4-2 would be used to link 
the user subroutine with the simulator. The subroutine would be on the sequential data referred 
to by the SYSIN DD statement of the FORTRANH procedure. 

Figures 4-4 through 4-12 are samples of the printed output created by the simulator. Figure 4-4a 
shows values of input parameters which vary with depth. In the following, NAMELIST names 
are given in parentheses. The first three columns are the layer index, the layer thickness (DZ), 
and the depth to the center of each layer. The latter is computed by the simulator using the input 
layer thicknesses. The next two columns are the initial wetness (WATERS) and the temperature 
(TEMPS) of each layer. The next three columns are the thermal conductivity of the solid material 
(TCONDS), porosity (PORSTY), and the volumetric heat capacity of the solids (VHCAPS). The 
next four columns are used in the model of the soil hydraulic properties, Equations 2-65 a and b. 
They are saturation wetness 0 $ (SATW), saturation matric potential \p s (SATP), saturation hydrau- 
lic conductivity K s (SATK), and soil texture parameter b (EB). The column labeled EM is the 
exponent 2b + 3 in the hydraulic conductivity model. The last column (ROOTS) is the root density 
profile, Pj from Equation 2-60. 
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** INPUT AND INITIALIZED PARAMETERS ** 


I THICKNESS 

DEPTH 

WETNESS 

TEMP 

TCONDS 


PORSTr 

VHCAP 

sat w 


S ATP 




satk 



bn 


2 M 



ROOTS 

10.500E 

00 

0.25 00 

00 

0.167E 

00 

0.2 79E 

03 

0. 250 E— 02 

0. 

A50E 

00 

0.500E 

00 

0.375E 

00 

- .142E 

02 

0 


122C-04 

C 


500E 

0 1 

0. 130F 

02 

C 

.0 

20.500E 

00 

0.7500 

00 

0. 193E 

00 

0 .260 E 

03 

0.250 E— 02 

0* 

45 CE 

00 

0 -500E 

00 

0.375E 

00 

- .142E 

02 

0 


12 2E-0 4 

0 


5C0E 

0 1 

0. 1 JOE 

02 

0 

.0 

30.200E 

0 1 

0.200D 

0 1 

0.229E 

00 

C.261E 

03 

C.250E-02 

0. 

A50E 

00 

0 .500E 

00 

0.375E 

00 

-.14 2E 

02 

0 


122b -04 

c 


5 OOE 

01 

0.133b 

02 

0 

.0 

AO* 2 00 E 

0 1 

0 . A 0 00 

0 1 

0 .2A3E 

00 

0 « 28 JE 

03 

0.250E-02 

0. 

45 0 E 

00 

0 .500E 

00 

0.375E 

00 

— .14 2E 

02 

c 


122E-04 

0 


5 COF 

01 

0.130b 

02 

0 

.0 

50 • 5 03 E 

01 

0.7500 

01 

0.256E 

00 

0.284E 

03 

0.250E-02 

0. 

450 E 

00 

0 .5005 

00 

0.375E 

00 

— . 14 2E 

02 

0 


12 2E-0 4 

0 


5C0E 

0 1 

0. 130E 

02 

0 

.0 

60.500E 

01 

0.1250 

02 

0.266E 

00 

0.285E 

03 

0.250E-02 

0. 

45 OE 

00 

0 ,500E 

00 

0.375E 

00 

— . 1 42E 

02 

0 


122E-04 

0 


SCOb 

0 1 

C. 1 JOE 

02 

0 

.0 

70* 5 00 E 

01 

0.1750 

02 

0 .266E 

00 

0-286E 

03 

0.250E-02 

0. 

45 OE 

CO 

0 « 500E 

00 

0.375E 

00 

- .142E 

02 

0 


12 2E-0 4 

0 


5 OOE 

01 

0. 130 b 

02 

0 

.0 

80. lOOE 

02 

0.2500 

02 

0.272E 

00 

0.286E 

03 

0.25 0 E— 02 

0. 

A50E 

00 

0 .500E 

00 

0.375E 

00 

— .14 2E 

02 

0 


12 2E-0 4 

0 


5 OOE 

01 

0. 1 30 E 

02 

0 

.0 

90. 1 00 E 

02 

0.35 00 

02 

0.265E 

00 

C.286E 

03 

0.250E-02 

0. 

45 CE 

00 

0 .500E 

00 

0.375E 

00 

-.14 2E 

02 

0 


122E-04 

0 


SCOE 

01 

0. 130E 

02 

0 

.0 

100. 100E 

02 

0.A50D 

02 

0.258E 

00 

0.286E 

03 

0. 250 E— 02 

0. 

45 OE 

00 

0.500E 

00 

0.375E 

00 

-.1 4 2E 

02 

0 


1 2 2b -0 4 

c 


500E 

01 

0. 130E 

02 

0 

.0 

l 10. 103E 

02 

0.5500 

02 

0 • 24 A E 

0 c 

C.287E 

03 

0.250E-02 

0. 

45 OE 

00 

0 .500E 

00 

0.375E 

00 

— .14 2E 

02 

0 


122E-04 

0 


SCOE 

01 

0. 130C 

02 

0 

.0 

120. lOOE 

02 

0.6500 

02 

0.239E 

00 

0.287E 

03 

0.250E-02 

0. 

45 OE 

0 0 

0 .5 OOE 

00 

0.375E 

00 

— .14 2E 

02 

c 


122E-04 

0 


500b 

01 

0. 13CE 

02 

c 

.c 

130. lOOE 

02 

0-75 00 

02 

0.21 AE 

oo 

0.287E 

03 

0.250 E— 02 

0. 

45 OE 

00 

0 .500E 

00 

0.375E 

00 

— .14 2E 

02 

0 


1 2 2b— C 4 

0 


5 COE 

01 

0. 1 30 E 

02 

0 

.0 

140. lOOE 

02 

0.8500 

02 

0. 170E 

00 

C.288E 

03 

0. 250 E— 0 2 

0. 

45 OE 

00 

0 .500E 

00 

0.375E 

00 

- .142E 

02 

c 


122E-04 

0 


5 0 0 E 

01 

0 a 1 30 E 

02 

0 

.0 

150. 1 OOE 

02 

0 .9500 

02 

0 • 153E 

00 

0.288E 

03 

0.250E-02 

0. 

450 E 

00 

0 .500E 

00 

0.375E 

00 

— • 1 4 2E 

02 

0 


122E-04 

c 


500b 

01 

C.130E 

02 

c 

.0 


Figure 4-4a. Input Parameters 
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Figure 4-4b is an example of the second output page, where the values of all other input parameters 
are given. The variable name appears first, followed by a brief description and its value. 

A typical page of table output is shown in Figure 4-5. The output time is in DDDHHMMSS.SS 
format, where DDD is the day (DDD = 0 is the first day), HH is the hour of the day, MM is minutes 
of the hour, and SS.SS is seconds. It is assumed that the simulation start time corresponds to mid- 
night. The first line also shows the current integrator step size in seconds. 

Next on this table are the values of some important variables related to moisture flow. Layer 
number, depth to the center of each layer, and wetness in the layer are shown in the first three 
columns. Next is the wetness flux in cm/sec (q 0 Equation 2-75a) at each of the NL + 1 layer 
boundaries. The first layer is the air-soil interface. The next two columns are the hydraulic con- 
ductivity (HYD COND) and matric potential, or pressure head (PHEAD) of the water in the layer. 
The next column lists values of the derivative of the matric potential with respect to depth 
Equation 2-73a). The I th derivative is evaluated at the boundary separating the I and I + 1 layers. 
The last entry in the column is not used, since the derivative is required only at the NL-1 interior 
boundaries. The column labeled DWDT is the derivative of wetness with respect to time , Equa- 
tion 2-54). If the sink term is zero, this is the spatial derivative of the wetness flux (— Equation 
2-7 2a). The last column is the value of the sink term Q, Equation 2-55. 

The next entry is the amount of water stored in the profile in centimeters. This will be calculated 
if the value of the input variable NWCUMS is greater than zero. 

The cumulative wetness flux at any of up to 10 layer boundaries can be computed. This is the 
integral of the flux over time, starting with the beginning of the simulation. The trapezoid rule is 
used to compute the integral, and the value of the integral is updated after each simulator time 
step. The number of integrated fluxes is set via NAMELIST parameter NWCUMS; the layer indices 
are input via the array IXWCUM. In the example of Figure 4-5 NWCUMS is 3, and the first three 
elements of IXWCUM are 1 , 2, and 16 (see Figure 4-4b). Boundary 1 is the air-soil interface, so 
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INTEGRATOR VARIABLES 
T STOP STOP TIME (SECONDS) 

WATERR ERROR TOLERANCE FOR WETNESS 
TEMERR ERROR TOLERANCE FOR TEMPERATURE 
ED ERROR WINDOW: WINDOW = (E/EDtE *ED) 

IFORCE FORCE MAX « STEP SIZE (O^NC. 1=YES) 

H INITIAL STEP SIZE (N/A IF IFORCE=t) 

HMAX MAXIMUM STEP SIZE (N/A IF IFQRCE=0) 


0.864000D 05 
0. 20E — 0 2 
C.20E 0 1 
0.20D 01 
1 

0.100000D 01 
0.1800000 04 


OUTPUT CONTROL VARIABLES 
D TOUT OUTPUT PERIOD (SECONDS) 

NFUNCT NUMBER OF GRAPHS OF WETNESS VS BEPTH 
AND TEMP VS DEPTH PER PLOT PAGE 
NWATRS NUMBER OF WETNESS VS TIME GRAPHS 
INOXW SOIL LAYERS FOR WETNESS PLOT 
NWFLUX NUMBER OF WETNESS FLUX VS TIME GRAPHS 
INDXWF BOUNDARIES FOR WETNESS FLUX PLOT 
N TEMPS NUMBER OF TEMP VS TIME GRAPHS 
INDXT SOIL LAYERS FORT TEMP PLOT 
NTFLUX NUMBER OF HEAT FLUX VS TIME GRAPHS 
INDXTF BOUNDARIES FORT HEAT FLUX PLOT 
IPR OUTPUT UNIT NUMBER FOR ALL PLOTS 

I TABLE PRINT TABLES < 0=NO, 1=YES) 

WL. WH FORCED MIN AND MAX FOR WETNESS PLOT 

TL* TH FORCED MIN AND MAX FOR TEMP* PLOT 

WFL.WFH FORCED MIN AND MAX FOR W FLUX PLOT 

TFL.TFH FORCED MIN AND MAX FOR H FLUX PLOT 

IF L AND H LIMITS ARE EQUAL * THEY ARE SET 


0.180000E 04 

6 

3 

1 230000000 

2 

^ 1 20000000 0 

^1 230 000 0 0 0 

1 250000000 

10 

1 

0*0 0*0 

0*0 0*0 

0*0 0*0 

0*0 0. 0 

THE PROGRAM 


NWCUMS NUMBER OF CUMULATIVE WETNESS FLUXES 3 

IXWCUM BOUNDARIES FOR CUMULATIVE WETNESS FLUXES 1 
NHCUMS NUMBER OF CUMULATIVE HEAT FLUXES 1 

IXHCUM BOUNDARIES FUR CUMULATIVE HEAT FLUXES l 

NDISK OUTPUT STATE TO DISK <0=NG*1=YES) 0 

IUDISK FORTRAN UNIT NUMBER FOR DISK OUTPUT 12 


OTHER MODEL PARAMETERS 

ITEMPS TEMPERATURE MODEL INDICATOR 

(0 = NO TEMP*, I = TEMP .PROFILE * 2= FORCE- PE ST ORB > 

I ROOTS INCLUDE ROOT MODEL (0=NO> 1=YES ) 

SPRES ROOT SPECIFIC RESISTANCE (SEC/C W) 

CPMIN PLANT CROWN POTENTIAL AT WILTING 

SFRAC FRACTION OF ET DEMAND TO SOIL E VAP 


0*1000000 07 
— * 1 50 00 OD 05 
0*1000000 00 


NEXT 3 PARAMETERS USED FOR GAUSSIAN TIME-DEPENDENT ET FUNCTI3-J WHEN ITEMPS = 0 

EMAX MAXIMUM RATE (CM/SEC) 0*300000E-04 

EDAY TOTAL DAILY E VAPOTRAN SP IR AT I ON (CM) 0.100000E 01 

EMAXT TIME OF MAX RATE (SEC SINCE START) 0*468000E 05 

THMIN LIMITING SURFACE WETNESS FGR E VAPOR ATI ON 0*5000 0 00-0 l 

RNTOT TOTAL RAIN FALL (CM) 0.0 

RNSTRT START TIME OF RAIN (SEC FRCM SIW START) 0*0 

RNSTOP STOP TIME OF RAIN (SEC FRCM SIM START)0.0 

RNRATE RAINFALL RATE (CM/SEC) 0*0 

JBOT BOTTOM WETNESS BOUND. COND. 1 

=0,FLUX-0.0; =1 .WETNESS IS CONSTANT; *2 «F LUX=CONDUCT I VI T Y 
THE FOLLOWING INPUT PARAMETERS ARE USED WHEN SOLVING THE HEAT BALANCE EQUATION (ITEMPS>0) 

ALB SURFACE SHORT WAVE ALBEDO 0.300000E 00 

CLOUDS FRACTIONAL CLOUD COVER 0.0 

CTRANS CLOUD TR A NSM I T T I V I TY TO SHCRT WAVES 0.200000E 00 

SUMJEC SUN DECLINATION ANGLE (DEGREES) 0.0 

XLAT LATITUDE 0.335000E 02 

TURB TURBIOITY FACTOR C.200000E 01 

CEVAPO COEFFICIENT FOR PET CALCULATION 0.0 

CEVAP1 COEFFICIENT FOR PET CALCULATION 0.500000E-06 

THESE ARE USED TO COMPUTE POTENTIAL E VAPO TRANSPIR AT 10 N AS A FUNCTION OF WIND SPEED (W, CM/SEC) 
AND VAPOR PRESSURE DIFFERENCE ( DVAP* MB.). THE MODEL IS PET = -(CEVAPO + CE V AP 1 *W * DVAP ) 


Figure 4-4b. 
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TIME (DDDHHMMSS) 


LUX 

I 

COMPUTED 

DEPTH 

AT 

1 

0.250000 

00 

2 

0 .750 COO 

00 

3 

0 .20000D 

01 

4 

0 .40 0000 

01 

5 

0.75 00 OD 

01 

6 

0 .12 50 00 

02 

7 

0 .17500D 

02 

a 

0 .250000 

02 

9 

0 .350000 

02 

10 

0.450000 

02 

1 1 

0.5500 00 

02 

12 

0.65000D 

02 

13 

0.750000 

02 

14 

0.3500 00 

02 

15 

0 .950000 

02 


1 20000.0 0 


WETNESS 

0*144490 CO 
0*194910 00 
0.22710D 00 
0*241150 00 
0.255320 00 
0.26476D 00 
0.267870 00 
0.27100D 00 
0.26613D 00 
0*2 586 50 00 
0.248360 00 
0.236380 00 
0.21449D 00 
0.177590 00 
0.153000 00 


FLUX 

-. 52689D— 05 
264520-05 
597360-06 
489690-06 
— . 36331 D— 06 
— • 23461 D-06 
0. 431550-00 
0. 63664D-07 
0. 26735D-06 
0. 26221 0-06 
0. 234970-06 
0. 18002D-06 
0. 191 500—06 
0. 172750-06 
0. 269710-07 
0. 269710-07 


HYD COND 

0.503270-10 
0.246460-03 
0.179800-07 
0.392380-07 
0.824270-07 
0. 1321 10-06 
0. 15380D-06 
0. 178830-06 
0. 14 1340-06 
0.975630-07 
0.575570-07 
0.30 2560-07 
0.855740-08 
0. 735310-09 
0. 1 0591D-09 


A TOR STEP 
IOUNOARY 
PHI AO 

SIZE (SECONDS) 
OPDZ 

= 0. 28125000000 C2 
DWDT SINK 

-.16/210 

04 

0.259560 

04 

- .524740-05 

0 .0 

- . 37 436 D 

03 

0 .1 60030 

03 

- .102390-05 

0 .0 

-.17 +320 

03 

0 *226 0 30 

02 

-.215330-07 

c .c 

-.12*120 

03 

0 .916220 

01 

-.252770-07 

0 .0 

— .9 7 0 54 D 

02 

0 .322090 

01 

— .12 870D— 07 

0 .0 

-.83 *500 

02 

0 .919250 

00 

-.238920-07 

0 .0 

-.76 J54D 

02 

0.573750 

00 

-.593480-00 

0 .0 

-.72 J51 D 

02 

-.682400 

00 

— *20 369 D— 07 

0 .c 

-•7Bi74D 

02 

-.120860 

01 

0 .51 388D — Ob 

0 .0 

-.93 *61 O 

02 

- .204700 

01 

0 .272390-07 

0 .0 

-•111 430 

03 

-.312690 

01 

0 .549530-07 

c .c 

-.14270 D 

03 

- .89239D 

01 

-.1 1482D-07 

0 .0 

-.231940 

03 

-.364170 

02 

0 .18751D-07 

0 .0 

-.59*11 D 

03 

-.659890 

02 

0 .145780-06 

0 .0 

-.125600 

04 

-.659890 

02 

0 .0 

0 .0 


CUMULATIVE WETNESS VARIABLES (CM CF WATER) 
TOT1L IN PROF ILE = O . 233 0 20 24 6 OD 02 


BOUNDARY 

1 

2 
1 6 


CUMULATIVE FLUX 
-.44 1422J7I4D-01 
-.33 19 12284 10-01 
0.8 3 5474 6 0 1 70-03 


EVAPOTRANSP IRATION OJTPUT - FLUXES 
TOTAL = - .526890-05 SOIL EVAP - 


IN CM/SEC 
-.526890-05 


PLANT TRA ISP IRATI UN = 0.0 


CROWN POTENTIAL (CM) = 0.0 


I 

D E 3 T H 


TEMP 


1 

0 .250000 

00 

0.30449D 

03 

2 

0 .75000D 

00 

0.300780 

03 

3 

0.200000 

01 

0.29525D 

03 

4 

0 .400000 

01 

0.289520 

03 

5 

0 .750000 

01 

0.284330 

03 

6 

0 .12 5000 

02 

0.283860 

03 

7 

0.17 5000 

02 

0.284720 

03 

a 

0 .250000 

02 

0.285630 

03 

9 

0.350000 

02 

0.286070 

03 

10 

0.4 5 00 OD 

02 

0.286450 

03 

1 1 

0.5 5 00 OD 

02 

0.286790 

03 

12 

0.550000 

02 

3.287120 

03 

13 

0.75 00 OD 

02 

0.28744D 

03 

1 4 

0.9500 00 

02 

0.2 87760 

03 

15 

0 .95 00 OD 

02 

0.288100 

03 


FLUX 

0. 74955D-02 
0. 10480D-01 
0. 638920-02 
0. 422150-02 
0. 21989D-02 
0. 140600-03 
2E890D-03 
185160-03 
661700-04 
57379D-04 
-. 51552D-04 
-. 480830-04 
-. 462510-04 
-. 454320-04 
451480-04 
-. 451 430-04 


TCOND 

0.13 8640- 02 
0. 14 3960-02 
0.146930-02 
0. 1 4778D-02 
0.149300-02 
0. 150930-02 
0.15 1620-02 
0.152330-02 
0. 15 152D-02 
0. 15 0250-02 
0. 148490-02 
O. 146450-02 
0. 142740-02 
0.136610-02 
0. 132670-02 


TDIrF 

0.981340-07 
C. 693570-07 
0.497970-07 
0.39761 0-07 
0.36a 780-07 
0 .403 100-07 
C .43o58D-07 
0.457820-07 
0.442320-07 
0.40 332 D— 07 
0.373 76 D— 07 
0.3* i 170-07 
0 .332400-07 
0.352 130-07 
C . 38 * 02 0—0 7 


WO IFF 

0.127420-1 1 
0 .862530-12 
0 .54546D-12 
0.360800-12 
0 .242420-12 
0. 223830-12 
0 .232500-12 
0 .242240-12 
0.255800-12 
0.272690-12 
0.293690-12 
0.317600-12 
0 .357230-12 
0. 421610-12 
0. 469210-12 


VH CAP 


0 .4 1 958 D 

00 

0 .46998D 

00 

0 .502170 

00 

0 .51621 0 

00 

0 .53 0380 

00 

0 .53981 0 

00 

0 .542920 

00 

0.54605D 

00 

0.54 1 190 

OC 

0 .53371 O 

00 

0 .523420 

00 

0 .51 1440 

OC 

0.48956D 

00 

0 .452670 

00 

0 .428C9D 

00 


OTDZ 

-.741840 01 
-.442000 01 
-.28648D 01 

-.148360 01 

- .9 36 620-01 
C.17114D 00 
0.121 9 3D OC 
0 .435540-01 
C .380280-01 
0 . J4513D-0 1 
0 .326060-01 
C .3 19880-01 
0 .325330-01 
0 .32S46D-0 1 
0 .335460-01 


DTEMPDT 

-. 1 42280- Cl 
0 .435250-02 
0. 86334D- C3 
0 • 78 3630- G3 
0.388090- 03 
0. 740070- 04 
-• 1 3 58 2D — 04 
-.21 790D-04 
-. 16244D-04 
-. 1 091 80-04 
— .662 7 JO - 05 
-.350260- C5 
-. 1 6734D-05 
-. 626590- 06 
0.0 


TERMS UF THE HEAT BALANCE EQUATION: FLUXES 
NET R A 01 A T I UN * O.I08938E-01 ET FLUX = 


IN CALv'CM**2/SEC 
— .308 758E— 0 2 


SINS IrJLE HEAT 


— .31 0 664 E— 0 3 


HEAT ABSORBED 


ITERATIVE SOLUTION OF HEAT BALANCE EQUATION GIVES SURFACE TEMPERATURE = 0.J0ba02O 03 OEGREES K 


BY SOIL 


0 .749553E-02 


CUMULATIVE HEAT VARIABLES (CAL/CM**2) 

TOTAL IN PROFILE = 0.14572405560 05 

BOUNDARY CUMULATIVE FLUX 

1 0.56 2062 761 4D 02 


Figure 4-5. Example of table output 



the first integrated flux is the integrated soil evaporation rate. The last entry is the amount of 
water lost through the bottom boundary. 

The evapotranspiration in cm/sec is shown in the next line. Values are given for the total flux 
(E, Equations 2-34 or 2-50), the soil evaporation flux (E s , Equation 2-63 if a root model is in- 
cluded; otherwise it equals the total flux), and the plant transpiration (E pl , Equation 2-64 if the 
root model is included). The crown potential (0 p , Equation 2-61 ) is also given here. 

The remaining output in this page gives values for parameters related to the temperature model. 

The first three columns are layer index, depth to center of the layer, and layer temperature (°K). 

The next column is the heat flux (q, Equation 2-75b) at layer boundaries. The next four columns 

are: thermal conductivity (X, Equation 2-23); thermal diffusivity (D T , Equation 2-20c); vapor 

wetness dif/usity (coefficient of ~~ in the definition of D , Equation 2-1 5c); and the volu- 

metric heat capacity (C, Equation 2-29). The derivative of temperature with respect to depth at 
dT 

the interior layers Equation 2-73 b) is given in the next column. Finally, the derivative of 

IT" 

temperature with respect to time ((-^p Equation 2-72b), which is used to integrate the state 
equations, is given in the last column. 

The values of the four terms of the heat balance equation (Equation 2-34) are given in the next 
line. These are net radiation (R, Equation 2-36), heat carried by the evapotranspiration flux (LE, 
Equation 2-42), sensible heat flux (H, Equation 2-46), and the heat absorbed by the soil (S, Equa- 
tion 2-34). The surface temperature (T s , Equation 2-47) is given in the next line. 

Finally, cumulative heat fluxes can be optionally computed. NHCUMS is the number of such 
fluxes, and the array IXHCUM contains the boundary indices. In this example, only the integrated 
surface flux is computed. 

The next six figures are examples of the plots generated by the simulator. The FORTRAN unit 
number of the output is the value of the NAMELIST variable IPR. 


4-16 



Figures 4-6 and 4-7 are moisture and temperature profiles respectively. To reduce the amount of 
printout, more than one profile can be put on a plot page. Variable NFUNCT is the number of such 
graphs per plot page; in this example NFUNCT = 6. If NFUNCT is zero, no profiles will be plotted. 
The correspondence between print character and output time is shown at the top of each page. If 
more than one print character fall in the same place, then the number of such characters is plotted. 
These two figures show that over the time interval from noon to 2:30 PM, the profiles change 
only in the surface layers. 

Figures 4-8 and 4-10 are plots of wetness and temperature respectively in selected layers as func- 
tions of time. NWATRS is the number of layers for which wetness is to be plotted, and INDXW 
is the array containing the indices of the chosen layers. In Figure 4-8 NWATRS is 3, and the in- 
dices are 1, 2, and 3. Similarly, NTEMPS is the number of layers for which temperature is to be 
plotted, and INDXT is the array containing the indices of the layers. In Figure 4-10 NTEMPS is 
3, and the first three elements of INDXT are 1 , 2, and 3. 

Figures 4-9 and 4-1 1 are plots of moisture and heat fluxes at selected boundaries. NWFLUX is 
the number of boundaries for which wetness flux is to be plotted, and the array INDXWF contains 
the indices of those boundaries. NTFLUX and INDXTF are the number of boundaries and their 
indices for the heat flux plot. 

The user can optionally control the vertical limits on all plots. NAMELIST variables WL and WH 
are respectively the lower and upper limits on wetness plots. (Figures 4-6 and 4-8). If WL = WH, 
then these values are ignored and appropriate limits are chosen from the data being plotted. Simi- 
larly, TL and TH are lower and upper limits for temperature plots (Figures 4-7 and 4-10); WFL and 
WFH are lower and upper limits for wetness flux plots (Figure 4-9); and TFL and TFH are lower 
and upper limits for heat flux plots. If any pair of lower and upper limits are equal, the actual 
limits are chosen from the data. 
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PLOT LABEL TO OUTPUT TIME CORRESPONDENCE (TIME IN OD OHHMM SS* i 

A 120 000 *00 E 1 23000.00 t 13 0300. CO C 

F 143 00 0.00 

*** WETNESS (CM3/CM3) VERSUS >EPTH (CM) *** 


1 33000. C 0 


0.2 72 
0.2 70 
0.266 
0.266 
0.264 
0.2 62 
0.260 
0.2 58 
0.2 56 
0.2 54 
0.2 52 
0.250 
0.2 40 
0.2 46 
0.244 
0.242 
0.240 
0.2 38 
0.2 36 
0.2 34 
0.232 
0.230 
.0.228 
0.226 
0.224 
0.2 22 
0.220 
0.2 18 
0.216 
0.2 14 
0.2 12 
0.210 
0.208 
0.206 
0.2 04 
0.202 
0.200 
0.193 
0. 1 95 
0. 1 94 
0. 1 92 
0. t 90 
0. 1 88 
0. 1 85 
0. 1 64 
0.182 
0. 1 60 
0. 1 78 
0. 1 76 
0. 1 74 
0. 1 72 
0. 1 70 
0.169 
0. 1 66 
0. 1 64 
0.162 
0. 160 
0. 1 58 
0. 1 56 
0.154 
0. 1 52 
0. 1 50 
0.148 
0.146 
0. 1 44 
0.142 
0.140 
C.l 33 
0.136 
0. 1 34 
0.132 
0. 1 33 
0 .1 28 
0. 1 25 
0. 1 24 
0.122 
0.120 
0.1 IS 










♦ 4 ♦ 4- * + 4 4. 


4 3 

♦ 3 


4 + — — 

48. 52. 


-- 4 - 

56. 


. 64 . 68. 


4. 4. 4-— — 

76. 80. 84. 


08. 


Figure 4-6. Example of plots of soil moisture profiles 
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PLOT L4 B EL tO OUTPUT TIME CORRESPCNDENCE (TIME 
A 120 000 *00 B 1 23000*00 

F 143000.00 


IN OODHHMMSS. J 

C 130000.00 


D 133000.00 


*** T EMPERAT JRE ( DEG K> VERSUS OEPTH (CM) *** 


E 14000C.GC 


307.2 4 4- 

306.8 2 

306.4 2 

306.3 F 

305.6 4 

305.2 +2 

304.8 4F 

304.4 A 

304.0 +C 

303.6 4 

303.2 4 

302.8 4 

302.4 ♦B 

302.0 4 

301.6 4 

301.2 ♦ F 

300.8 +A E 

300.4 4 

300.0 4 

299.6 4 D 

299.2 4 

299.8 4 

298.4 4 C 

298.0 4 

297.5 4 

297.2 4 

296.8 4 B 

296.4 4 

296.0 4 

295.5 4 F 

295.2 4 A 

294.8 4 

294.4 4 E 

294.0 4 

293.6 4 D 

293.2 4 

292.9 4 

292.4 4 C 

292.0 4 

291.6 4 

291.2 4 

290.8 4 B 

290.4 4 

290.0 4 

289.6 4 A 

289.2 4 

288.8 4 

288.4 4 

288.0 4 

287.6 4 

287.2 4 

286.9 4 

286.4 4 

286.0 4 

285.6 4 

285.2 4 

284.8 4 

284.4 4 
2e4.0 4 

283.6 4 4 - 

0. 4 



Figure 4-7. Example of plots of soil temperature profiles 
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*** * ET NESS ViRSUS TIME * ** 


START TIME (DDDHHMMSS.) = 0*0 

PRINT CHARACTER TO SOIL LEVEL CO RflESPO hDENCE 
At 9 2 C 3 

0*2295 C — C + ♦ + 4 « ♦ - 

0*2200 «■ C CC CC CC CC CC CC CC CC 

0*2265 ♦ 

0.2250 ♦ 

0.2235 ♦ 

0*2220 ♦ 

0.2205 ♦ 

0.2190 ♦ 

0.2175 ♦ 

0.2160 ♦ 

0.2145 4 
0. 2 1 30 + 

0.2115 4 
0.2100 «- 
0.2085 «• 

0.20 70 ♦ 

0.2055 ♦ 

0.2040 ♦ 

0.20 25 ♦ 

0.20 10 «• 

0.1995 ♦ 

0.1980 ♦ B 

0.1965 ♦ B B E B 

0.1950 ♦ B 8 B B 

0.1935 B B a B b B 
0.1920 ♦ 

0.19 05 ♦ 

0.1890 ♦ 

0.1875 «■ A 

0.1860 ♦ AAA 

0.1845 4 - A A 

0.1830 ♦ A 

0.1015 ♦ A 

0.1800 ♦ A 

0.1705 ♦ A 

0.1770 ♦ A 

0.1755 ♦ A 

0.1740 ♦ 

0.1725 ♦ A 
0.1710 ♦ 

0.1695 ♦ A 
0.16 00 «■ 

0.1665 A 
0.1650 ♦ 

0.1635 ♦ 

0.1620 ♦ 

0.1605 -#• 

0.1590 ♦ 

0.15 75 ♦ 

0.15 60 
0.1545 4 
0.1530 4 

0.1515 ♦ 

0.1500 ♦ 

0. 1485 ♦ 

0.14 70 ♦ 

0.1455 ♦ 

0.1440 ♦ 

0.1425 4 

0.1410 ♦ 

0.1395 ♦ 

0.1300 + 

0.1365 + 

0.1350 
0.1335 4 
0.1320 + 

0.1305 ♦ 

0.1290 ♦ 

0.12 75 + 

0.1260 
0.1245 ♦ 

0.1230 + 

0.1215 «• 

0.1200 + 

0.1185 + 

0.1170 4 ♦ + *• + ♦ * *■- 


TIME AXIS IS IN HOURS SINCE START 


CC CC CC CC CC CC CC 


C C C C C C 


b a a a b b 


a as as bb b b 
a a E B B 

B 8 B B B 


A A 

AAA A 


Figure 4-8. Example of plots of soil moisture versus time in selected layers 
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*** WETNESS FuUX (CM/SEC) VERSJS TIME *** 

START TIME (DDDHHMMSS.) = 0.0 TIME AXIS IS IN HOURS SINCE START 

PRINT CHARACTER TO SOIL LEVEL CORRESPONDENCE 
A 1 B 2 

Y-AXIS VALUES MUST BE MULTIPLIED BY l*0E-05 TO OBTAIN CORRECT VALUE 


■r 

to 


-O. 0 03 
-0.0 16 
-0.0 24 
-0 . 0 32 
-0 * 0 40 
-0.040 
-0.0 56 
-0.0 64 
-0.0 72 
-0.080 
-0.0 80 
-0.0 96 
-0.104 
- 0.112 
-O. 1 20 
-O. 1 28 
-O. 1 36 
-0. 1 44 
-0. 1 52 
-0.160 
-0.168 
-O. 1 76 
-0 . 1 84 
-0 . 1 92 
-0.200 
-0 .2 03 
-0.216 
-0.224 
-0.2 32 
-0.240 
-0.248 
-0.256 
-0.2 64 
-0 . 2 72 
-0.2 80 
-O .288 
-0.2 96 
-0.304 
-0.312 
-O. 320 
-0 .329 
-0 .3 36 
-0.344 
-0 • 3 52 
-0.3 60 
-0.368 
-0 . 3 76 
-0.384 
-0.3 92 
-0.4 00 
-0.408 
-0.4 16 
-0.424 
-0 . 4 32 
-0.440 
-0.448 
-0 . 4 56 
-O . 4 64 
-0.4 72 
-0.460 
-0.4 08 
-0.4 96 
-0. 504 
-0.5 12 
-0.520 
-0.523 



Figure 4-9. Example of plot of soil moisture flux versus time at selected boundaries 
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*** TEMPERAT J IL ( DEG K) VERSUS TIME *** 


START TIME (DDDHHMMSS.) = 0.0 

PRINT CHARACTER TU SOIL LEVEL CORRESPONDENCE 
A l B 2 C 3 

307 .3 f- — — — 4 .— — — f- - — — — 4—— - — 4—— — — 4— — — — 4— — — — 4 — — — — - 

306.5 4 

306.0 4 

305.5 4 

305.0 4 

304.5 4 

304.0 4 

303.5 4 

303.0 4 

302.5 4 

302.0 4 
301 .5 4 

301.0 4 

300.5 4 

300.0 4 

299.5 4 

299.0 4 

298.5 4 

298.0 4 

297.5 4 

297.0 4 

296.5 4 

296.0 4 

295.5 4 A | 

295.0 4 

294.5 4 

294.0 4 

293.5 4 

293.0 4 

292.5 4 n 

292.0 4 A U 

291.5 4 

291 .0 4 

290 .5 4 

290.0 4 

289.5 4 ( 

289.0 4 B 

208 .5 4 

288.0 4 A 

287.5 4 

287.0 4 c 

286.5 4 R 

286.0 4 

285.5 4 

285.0 4 c 

284.5 4 A 

284.0 4 

203.5 4 Q 

283.0 4 

282.5 4 r 

282.0 4 A 

281.5 4 

28 1.0 C C 

280.5 4 C B C 

280.0 H C 

279.5 4 C C 

279.0 A 2 D C C 

278.5 

278.0 

277.5 

277.0 

276.5 

276.0 

275.5 

275.0 

274.5 
274 .0 


TIME AXIS 15 IN HOURS SINCE START 


-4- — -4- -A-A----4 4-- — 4- ♦ 4 — f f + - -4 

A A A 

A 

B B 

8 

A A 

8 8 

B 


□ 


A B 

A 


C C 


c c 


c c 


c c 

A 

B B 




2 . 


7. 


10 . 


12 . 


14. 


17. 


- 4 4 4 4 4- 

19. 20. 21. 22. 23. 


4 

24. 


Figure 4-10. Example of plots of temperature versus time in selected layers 


< 
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level correspondence 

Y-AXIS VALUES MUST BE MULTIPLIED BY I.OE-Ol TO OBTAIN CORRECT VALUE 


«** HEAT - LU< (CAL/CM**2/SEC> VERSJS TIME 
TIME AXIS IS in HOURS SINCE START 


0 . 1 03 
0. 1 04 
0*100 
0*096 
0.0 92 
0*0 89 
0.0 64 
0. 080 
0. 0 76 
0.0 72 
0 .068 
0.0 64- 
0. 060 
0.0 56 
0.0 52 
0.048 
0.044 
0.040 
0.036 
0.032 
0.0 28 
0.0 24 
0.0 20 
0.016 
0.012 
0.008 
0.0 04 
0.0 

-0.004 
-0.008 
- 0.012 
-0.0 16 
- 0.0 20 
-0.0 24 
-0.028 
-0.0 32 
-0.0 36 
-0.040 
-0.0 44 
-0.0 43 
-0.052 
-0 . 0 55 


V 4 


♦ 

4 

4 

4 

4 

4 

♦ B 

4 

4 

c c c 

4 

B 

AAA 

4 

4 

♦ B 
4 

+ 

4 

4 

4 

0. 1« 






-4 4 - 4 - -- 4 - ♦“ 

B B 


«. 4 4 4 - - 4 -- 4 


4 4 4 4 


B 

/k A A 


A A A 


A A 
B 


C C 


A 

B 

C C 

A 


B 

c c c 


N A c 

B C C 

CCCCCCC2CC 


C C 

B C 


8 0 

A 


2 C 


C c C c 

B 


AA A 2 2 A AA A A 

e e 

B 0 

B 

B 




e 

BA A 
A 


2 A 2 

A A B 


C C 


2 2 


15. 


+ 4 

16. 17. 


19. 


20 . 


21. 22. 23. 


Figure 4-1 1 . Example of plots of heat flux versus time at selected layer boundaries 



Moisture and temperature profiles can also be output to disk or tape. The NAMELIST variable 
NDISK controls whether this is done, and IUDISK is the FORTRAN unit number of the DD card 
which points to the output data set. This can be a sequential data set on disk or tape, or a member 
of a partitioned disk data set. 

The output records are unformatted. The record length is 8* (NL + 1) bytes, where NL is the num- 
ber of soil layers. This is 2* NL + 2 words of data per record. The records are written with sub- 
routine FWRITE of the FTIO package. They can be read with subroutine FREAD from the same 
package. If desired, a FORTRAN read statement with a character format can be used instead. 

One header record and one data record at each simulator output time are created. The first record 
contains the number of soil layers and thickness of each layer. Each data record contains the out- 
put time, the wetness profile, and the moisture profile. Record formats are shown in Figure 4-12. 
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Figure 4-12. Record format for profile data set created by the simulator 
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